1 Sec. 2.2
- 1.1 Box-Ljung test (IBM and CRSP Value-Weighted Index)
- 1.2 Find IBM monthly and yearly expected return
2 Sec. 2.4
3 Sec. 2.5
4 Sec. 2.6
- 4.1 3M log price (Sec. 2.6.3)
- 4.2 Find 3M monthly and yearly expected return
5 Sec. 2.7
6 Sec. 2.8
- 6.1 J&J quarterly earnings
- 6.2 CRSP Decile 1 Index
7 Sec. 2.9
- 7.1 Weekly 1yr and 3yr interest rate
  - 7.1.1 Exploration
  - 7.1.2 Fit models
8 Sec. 2.10
- 8.1 Example of consistent estimation for covariance matrix (1yr and 3yr rates)
- 8.2 Prove that the equation on pp101 is the $j^{th}$ diagnoal element of the equation (2.49).
9 Sec. 2.11
- 9.1 Daily simple return of value-&equal-weighted index from CRSP

# install.packages("bizdays")
# install.packages("RQuantLib")
# install.packages("fUnitRoots")
# install.packages("TSA")
# install.packages("devtools")
# install.packages("Dict")

# install.packages("devtools")

also installing the dependencies ‘askpass’, ‘credentials’, ‘openssl’, ‘sys’, ‘zip’, ‘gitcreds’, ‘httr2’, ‘ini’, ‘httpuv’, ‘mime’, ‘xtable’, ‘fontawesome’, ‘sourcetools’, ‘later’, ‘promises’, ‘jquerylib’, ‘sass’, ‘systemfonts’, ‘textshaping’, ‘tinytex’, ‘diffobj’, ‘rematch2’, ‘clipr’, ‘gert’, ‘gh’, ‘purrr’, ‘rappdirs’, ‘rprojroot’, ‘rstudioapi’, ‘whisker’, ‘cachem’, ‘shiny’, ‘callr’, ‘processx’, ‘bslib’, ‘downlit’, ‘httr’, ‘ragg’, ‘rmarkdown’, ‘xml2’, ‘htmlwidgets’, ‘stringr’, ‘prettyunits’, ‘xopen’, ‘brew’, ‘commonmark’, ‘stringi’, ‘cpp11’, ‘brio’, ‘praise’, ‘ps’, ‘waldo’, ‘usethis’, ‘desc’, ‘fs’, ‘memoise’, ‘miniUI’, ‘pkgbuild’, ‘pkgdown’, ‘pkgload’, ‘profvis’, ‘rcmdcheck’, ‘remotes’, ‘roxygen2’, ‘rversions’, ‘sessioninfo’, ‘testthat’, ‘urlchecker’


Warning message in install.packages("devtools"):
“installation of package ‘gert’ had non-zero exit status”
Warning message in install.packages("devtools"):
“installation of package ‘usethis’ had non-zero exit status”
Warning message in install.packages("devtools"):
“installation of package ‘devtools’ had non-zero exit status”
Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

# install.packages("Dict")

also installing the dependencies ‘tidyselect’, ‘dplyr’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

# devtools::create("AFTSCode")

library(fUnitRoots)
library(bizdays, RQuantLib)
library(TSA) # eacf
library(forecast)
print(packageVersion("bizdays"))
print(packageVersion("RQuantLib"))
library(knitr)
library(IRdisplay)
library(fBasics)
library(Dict)
library(sandwich) # Use heteroskedasticity-consistent (HC) variance estimation
library(lmtest) # Apply HC variance to linear models
library(tseries) # Jarque-Bera test

Attaching package: ‘bizdays’


The following object is masked from ‘package:stats’:

    offset


Attaching package: ‘TSA’


The following objects are masked from ‘package:stats’:

    acf, arima


The following object is masked from ‘package:utils’:

    tar


Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Registered S3 methods overwritten by 'forecast':
  method       from
  fitted.Arima TSA 
  plot.Arima   TSA 


Attaching package: ‘forecast’


The following object is masked from ‘package:bizdays’:

    bizdays

[1] ‘1.0.16’
[1] ‘0.4.21’

Attaching package: ‘fBasics’


The following objects are masked from ‘package:TSA’:

    kurtosis, skewness


Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

# detach("package:AFTSCode", unload = TRUE) # First detach the package
# unloadNamespace("AFTSCode")

library(AFTSCode)

# find("plot_pacf_acf"); find("draw_arima_forecast_fig")

# install_redir_output <- function(lib, log_fname) {
#     # Open a file connection for writing
#     sink(log_fname, append = TRUE)

#     # Install packages (output will be redirected to the file)
#     install.packages(lib)

#     # Close the file connection
#     sink()
# }

# plot_time_fig <- function(ts, main, xlab, ylab) {
#     par(bg = "white")
#     plot(ts, type = 'l', main = main, xlab = xlab, ylab = ylab)
#     points(ts, pch = '*')
# }

# cal_phi_0 <- function(arima_md, ord, digits=6) {
#     phi_0 = as.numeric((1-sum(arima_md$coef[1:ord]))*arima_md$coef['intercept'])
#     print(format(phi_0, digits = digits))
#     phi_0
# }

# get_mu <- function(arima_mod, digits=6) {
#     mu = as.numeric(arima_mod$coef['intercept'])
#     print(format(mu, digits = digits))
#     mu
# }

# plot_arima_forecast_fig <- function(
#     da_ts, eotr, h, npts, frequency, 
#     order, fixed, method, transform.pars,
#     main, xlab, ylab, ylim = NULL
# ) {
#     par(bg = "white")
#     # arima model
#     tr_da_ts = ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
#     if (is.null(transform.pars)) {
#         ts_fm3 = arima(tr_da_ts, order = order, fixed = fixed, method = method)
#     } else {
#         ts_fm3 = arima(tr_da_ts, order = order, fixed = fixed, method = method, transform.pars = transform.pars)
#     }
#     print(ts_fm3$nobs)
#     # Forecast
#     ts_fm3$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
#     ts_fc_res = forecast(ts_fm3, h = h)
#     # Plot forecast
#     xmin = time(da_ts)[eotr]-npts/frequency
# #     xmax = end(da_ts)[1]+(end(da_ts)[2]+2)/frequency
#     xmax = time(da_ts)[eotr]+(max(h, length(da_ts)-eotr)+1)/frequency
#     cat(xmin, ";", xmax)
    
#     plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
#     # Plot forecast mean
#     dummy_1st_fmean_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_fmean_ts)
#     points(dummy_1st_fmean_ts, pch = 1)
#     # Plot confidence interval
#     dummy_1st_flower_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[,2])), frequency = frequency, start = end(tr_da_ts))
#     dummy_1st_fupper_ts = ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[,2])), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_flower_ts, lty=2)
#     lines(dummy_1st_fupper_ts, lty=2)
#     # Plot original data
#     orig_plot_ts = ts(da_ts[(eotr-npts+1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr]-(npts-1)/frequency)
#     lines(orig_plot_ts)
#     points(orig_plot_ts, pch = 19)
#     ts_fc_res
# }

# comb_forecast_res <- function(forecast_obj, da_ts, eotr, freq) {
#     display(summary(forecast_obj))
#     fc_std_err = (forecast_obj$upper[,2]-forecast_obj$lower[,2])/2/qnorm(p = 0.975)
#     actual_ts = ts(da_ts[(eotr+1):length(da_ts)], frequency = freq, start = time(da_ts)[eotr+1])
#     display(forecast_obj$mean); display(fc_std_err); display(actual_ts)
#     multistep_ahead_forecast_tb = cbind(forecast_obj$mean, fc_std_err, actual_ts)
#     dimnames(multistep_ahead_forecast_tb)[[2]] <- c("Forecast", "Std. Error", "Actual")
#     multistep_ahead_forecast_tb
# }

# plot_unit_root_figs <- function(da_ts, freq, start, lag_max, main, xlab, ylab) {
#     par(bg = "white")
# #     lgdp_ts = ts(lgdp, frequency = 4, start = c(1947 ,1))
#     diff_da_ts = ts(diff(da_ts), frequency = freq, start = start)
#     plot(da_ts, col = 'black', type = 'l', main = main, xlab = xlab, ylab = ylab)
#     acf(da_ts, lag.max = 20)
#     plot_time_fig(diff_da_ts, main = main, xlab = xlab, ylab = paste("diff(", ylab, ")", sep = ""))
#     pacf(diff_da_ts, lag.max = 20)
# }

# plot_acf <- function(da, lag.max = NULL, main=NULL, w=NULL, h=NULL, ...) {
#     if (is.null(w) | is.null(h)) {
#         par(bg = 'white')
#     } else {
#         par(bg = 'white', pin = c(w, h))
#     }
#     plot(acf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
# }

# plot_pacf_acf <- function(da, lag.max = NULL, main=NULL, w=NULL, h=NULL, ...) {
#     if (is.null(w) | is.null(h)) {
#         par(mfrow = c(2, 1), bg = 'white')
#     } else {
#         par(mfrow = c(2, 1), bg = 'white', pin = c(w, h))
#     }
#     plot(pacf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
#     plot(acf(da, lag.max = lag.max, plot = F, `drop.lag.0` = F, ...), main = main)
# }

Sec. 2.2¶

Box-Ljung test (IBM and CRSP Value-Weighted Index)¶

IBM monthly returns shows almost no serial correlation.
CRSP Value-Weighted Index monthly returns shows medium momentum at lag-1 and medium reverting at lag-3.

da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header=T)
da[1:5,]

find("acf")

# IBM (Fig. 2.1)
plot_pacf_acf(da$ibmrtn, lag.max = 100)

# CRSP Value-Weighted Index (Fig. 2.2)
plot_pacf_acf(da$vwrtn, lag.max = 100)

Box.test(da$ibmrtn, lag = 5, type = 'Ljung')

	Box-Ljung test

data:  da$ibmrtn
X-squared = 3.3682, df = 5, p-value = 0.6434

libm = log(da$ibmrtn+1)
Box.test(libm, lag = 5, type = 'Ljung')

	Box-Ljung test

data:  libm
X-squared = 3.5236, df = 5, p-value = 0.6198

Find IBM monthly and yearly expected return¶

For the log(Return)
- $\mu_{monthly}\approx 0.010891$
- $\sigma_{monthly}\approx 0.070333$

# µ, sigma, and monthly price ratio
ibm_basic_stats <- apply(as.data.frame(libm), 2, basicStats)
ibm_basic_stats$libm['Mean',];
ibm_basic_stats$libm['Stdev',];
mnthly_ibm_pr_ratio = exp(ibm_basic_stats$libm['Mean',]+ibm_basic_stats$libm['Stdev',]^2/2)
mnthly_ibm_pr_ratio; mnthly_ibm_pr_ratio^12;

Sec. 2.4¶

AR(3)¶

gnp = scan(file = "../AFTS_data/Ch02/dgnp82.txt")
gnp1 = ts(gnp, frequency = 4, start = c(1947, 2))

plot_time_fig(gnp1, main = "Gross National Product", xlab = "Year", ylab = "Growth")

Find the AR order using AIC

m1 = ar(gnp, method = 'mle')
m1$order

m2 = arima(gnp, order = c(3, 0, 0))
m2

Call:
arima(x = gnp, order = c(3, 0, 0))

Coefficients:
         ar1     ar2      ar3  intercept
      0.3480  0.1793  -0.1423     0.0077
s.e.  0.0745  0.0778   0.0745     0.0012

sigma^2 estimated as 9.427e-05:  log likelihood = 565.84,  aic = -1121.68

$\phi_0 = \mu*(1-\phi_1-\phi_2-\phi_3)$

cal_phi_0(m2, ord=3)

[1] "0.00472311"

Residual standar error

sqrt(m2$sigma2)

Characteristic equation

?polyroot

p1 = c(1, -m2$coef[1:3])
roots = polyroot(p1)
roots

for (i in 1:3) {
    print(i)
    print(as.numeric(Mod(1-m2$coef[1]*roots[i]-m2$coef[2]*roots[i]^2-m2$coef[3]*roots[i]^3)))
}

[1] 1
[1] 2.237726e-16
[1] 2
[1] 4.641385e-16
[1] 3
[1] 2.237726e-16

Modulus of the roots

roots; Mod(roots)

1/roots; Mod(1/roots)

Compute the average length of business cycles
Doesn't matter if using the characteristic roots or the inverse of them.
Unit is "quarter". So average business cycle is about 10.66 quarters.

k = 2*pi/acos(0.434406493995099/0.522654950320849)
k

k = 2*pi/acos(1.59025281352281/1.91330819575347)
k

AIC¶

gnp = scan(file = '../AFTS_data/Ch02/dgnp82.txt')
ord = ar(gnp, method = 'mle')
format(ord$aic, digits = 6)

ord$order; ord$order.max

AR(3) example¶

da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]

Parameter estimation for value-weighted index¶

Try simple return

vwrtn_m3 = arima(da$vwrtn, order = c(3, 0, 0))
vwrtn_m3

Call:
arima(x = da$vwrtn, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1158  -0.0187  -0.1042     0.0089
s.e.  0.0315   0.0317   0.0317     0.0017

sigma^2 estimated as 0.002875:  log likelihood = 1500.86,  aic = -2991.73

vwrtn_ts = ts(da$vwrtn, frequency = 12, start = c(1926, 1))
plot_time_fig(vwrtn_ts, main = "Value-weighted index return", xlab = "Year", ylab = "Return")

est_mrtn = get_mu(vwrtn_m3)
(1+est_mrtn)^12-1

[1] "0.00894788"

Try log-return

lvwrtn_ts = ts(log(1+da$vwrtn), frequency = 12, start = c(1926, 1))
lvwrtn_m3 = arima(lvwrtn_ts, order = c(3, 0, 0))
lvwrtn_m3

Call:
arima(x = lvwrtn_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1105  -0.0139  -0.0953     0.0074
s.e.  0.0315   0.0317   0.0318     0.0017

sigma^2 estimated as 0.002886:  log likelihood = 1498.89,  aic = -2989.78

plot_time_fig(lvwrtn_ts, main = "Value-weighted index log-return", xlab = "Year", ylab = "log(Return)")

est_mlrtn = get_mu(lvwrtn_m3)
exp(est_mlrtn*12)-1

[1] "0.00744436"

Calculate the actual yearly return¶

The 9.3% per annum is close to 11.3% per annum from the estimated model on simple return.
If I transfer the simple return to log return using log(1+rtn), the estimated return per annum is 9.3%, which is closer to the actual value.

vw = da$vwrtn # Value-weighted index
t1 = prod(vw+1) # Simple gross return
t1; length(vw); t1^(12/length(vw))-1; # monthly return per annum

Model checking¶

m3 = arima(da$vwrtn, order = c(3, 0, 0))
m3

Call:
arima(x = da$vwrtn, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1158  -0.0187  -0.1042     0.0089
s.e.  0.0315   0.0317   0.0317     0.0017

sigma^2 estimated as 0.002875:  log likelihood = 1500.86,  aic = -2991.73

Compute the intercept $\phi_0 = \mu*(1-\phi_1-\phi_2-\phi_3)$

cal_phi_0(m3, ord=3)

[1] "0.00901194"

Compute the standard error of residuals

attributes(m3)

c(
    sqrt(sum(m3$residuals^2)/(m3$nobs-2*length(m3$coef)-1)),
    sqrt(sum(m3$residuals^2)/(m3$nobs-length(m3$coef))),
    sqrt(sum(m3$residuals^2)/m3$nobs)
)

sqrt(m3$sigma2)

Ljung-Box test (R uses 12 dof)

Box.test(m3$residuals, lag = 12, type = 'Ljung')

	Box-Ljung test

data:  m3$residuals
X-squared = 16.352, df = 12, p-value = 0.1756

Compute p-value using 9 dof (12-3=9, the parameter for AR does not include the intercept.)

pv = 1-pchisq(16.35,9)
pv

Refine the model to rule out statistically insignificant coef by fixing the AR(2) coef to zero (pp51)

mm3 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA))
mm3

Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)):
“some AR parameters were fixed: setting transform.pars = FALSE”

Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA))

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38

mm3_1 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "ML")
mm3_1; cal_phi_0(mm3_1, ord = 3)

Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”

Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "ML")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38

[1] "0.0088825"

mm3_2 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "CSS")
mm3_2; cal_phi_0(mm3_2, ord = 3)

Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”

Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "CSS")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1126    0  -0.1064     0.0091
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002878:  part log likelihood = 1500.33

[1] "0.00899938"

mm3_3 = arima(da$vwrtn, order = c(3,0,0), fixed = c(NA,0,NA,NA), method = "CSS-ML")
mm3_3; cal_phi_0(mm3_3, ord = 3)

Warning message in arima(da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), :
“some AR parameters were fixed: setting transform.pars = FALSE”

Call:
arima(x = da$vwrtn, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA), method = "CSS-ML")

Coefficients:
         ar1  ar2      ar3  intercept
      0.1136    0  -0.1063     0.0089
s.e.  0.0313    0   0.0315     0.0017

sigma^2 estimated as 0.002876:  log likelihood = 1500.69,  aic = -2993.38

[1] "0.00888153"

cal_phi_0(mm3, ord = 3)

[1] "0.00888153"

In the arima software, when estimating the residual variance, they use sample size $T$, instead of $T-2*p-1$.

sum(mm3$residuals^2)/(length(da$vwrtn)); mm3$sigma2; sqrt(mm3$sigma2)

box_test <- Box.test(mm3$residuals, lag = 12, type = 'Ljung')
box_test

	Box-Ljung test

data:  mm3$residuals
X-squared = 16.828, df = 12, p-value = 0.1562

attributes(box_test)

Compute p-value using 10 dof (12-2=10, the parameter for AR does not include the intercept.)

pv = as.numeric(1-pchisq(box_test$statistic, 10))
pv; box_test$statistic

AR(3) forecasting example¶

See issues related to scoping for forecast: link

# install.packages("forecast")

also installing the dependencies ‘farver’, ‘labeling’, ‘munsell’, ‘R6’, ‘RColorBrewer’, ‘viridisLite’, ‘pkgconfig’, ‘xts’, ‘TTR’, ‘curl’, ‘gtable’, ‘isoband’, ‘scales’, ‘tibble’, ‘withr’, ‘quadprog’, ‘quantmod’, ‘colorspace’, ‘fracdiff’, ‘generics’, ‘ggplot2’, ‘lmtest’, ‘magrittr’, ‘tseries’, ‘RcppArmadillo’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]

f_da_ts = ts(da$vwrtn[1:tr_n], frequency = 12, start = c(1926, 1))
ts_fm3 = arima(f_da_ts, order = c(3,0,0))
# ts_fc_res = forecast(ts_fm3, h = 12)
ts_fc_res = predict(ts_fm3, n.ahead=12)
ts_fc_res

f_da_ts = ts(da$vwrtn[1:tr_n], frequency = 12, start = c(1926, 1))
ts_fm3 = arima(f_da_ts, order = c(3,0,0))
ts_fm3$x <- f_da_ts
ts_fc_res = forecast(ts_fm3, h = 12)
summary(ts_fc_res)

Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = f_da_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2969.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715

cal_phi_0(ts_fm3, ord = 3)

[1] "0.00972842"

sqrt(ts_fm3$sigma2)

# Use forecast
par(bg = "white")
xend = end(f_da_ts)[1]+end(f_da_ts)[2]/12
xmin = xend-1
xmax = xend+1
plot(ts_fc_res, xlim = c(xmin, xmax))

# Use predict
par(bg = "white")
xend = end(f_da_ts)[1]+end(f_da_ts)[2]/12
xmin = xend-1
xmax = xend+1
plot(ts_fc_res$pred, xlim = c(xmin, xmax), ylim = c(-0.3, 0.4))

attributes(ts_fc_res)

Make multi-steps forecast plots (Fig. 2.7)

freq = 12
vwrtn_ts = ts(da$vwrtn, frequency = freq, start = c(1926, 1))
eotr=984
ts_fc_res = plot_arima_forecast_fig(
    da_ts=vwrtn_ts, 
    eotr=eotr, h=12, npts=12, frequency=freq, 
    order=c(3,0,0), fixed=NULL, method=NULL, transform.pars=TRUE,
    main="Forecasts from ARIMA(3,0,0) with non-zero mean for ewrtn (Fig. 2.7)", 
    xlab="Year", ylab="Return", ylim=c(-0.21, 0.21)
)

2006.917 ; 2009

Forecast values (Table 2.2)

summary(ts_fc_res)

Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = c(3, 0, 0))

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2967.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715

install_redir_output(lib = "knitr", log_fname = "knitr_install.log")

also installing the dependencies ‘highr’, ‘xfun’, ‘yaml’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

vwrtn_multistep_ahead_forecast_tb <- comb_forecast_res(ts_fc_res, vwrtn_ts, eotr, freq)
vwrtn_multistep_ahead_forecast_tb

Forecast method: ARIMA(3,0,0) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, fixed = fixed, method = method)

Coefficients:
         ar1      ar2      ar3  intercept
      0.1034  -0.0201  -0.1089     0.0095
s.e.  0.0317   0.0318   0.0317     0.0017

sigma^2 estimated as 0.00284:  log likelihood = 1488.85,  aic = -2967.71

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.234691e-05 0.05328902 0.03757996 101.6395 166.5597 0.7001715
                    ACF1
Training set 0.002856914

Forecasts:
         Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
Jan 2008    0.007455686 -0.06083694 0.07574831 -0.09698887 0.1119002
Feb 2008    0.015952488 -0.05270459 0.08460957 -0.08904946 0.1209544
Mar 2008    0.011701663 -0.05695841 0.08036174 -0.09330486 0.1167082
Apr 2008    0.009806653 -0.05927763 0.07889093 -0.09584863 0.1154619
May 2008    0.008770886 -0.06033071 0.07787249 -0.09691089 0.1144527
Jun 2008    0.009164644 -0.05993698 0.07826627 -0.09651718 0.1148465
Jul 2008    0.009432520 -0.05967458 0.07853962 -0.09625768 0.1151227
Aug 2008    0.009565095 -0.05954249 0.07867268 -0.09612583 0.1152560
Sep 2008    0.009530555 -0.05957703 0.07863814 -0.09616037 0.1152215
Oct 2008    0.009495151 -0.05961250 0.07860281 -0.09619589 0.1151862
Nov 2008    0.009477748 -0.05962992 0.07858541 -0.09621331 0.1151688
Dec 2008    0.009480420 -0.05962725 0.07858808 -0.09621063 0.1151715

class(vwrtn_multistep_ahead_forecast_tb); attributes(vwrtn_multistep_ahead_forecast_tb)

Find Volume-weighted monthly and yearly expected return¶

For the log(Return)
- $\mu_{monthly}\approx 0.007403$
- $\sigma_{monthly}\approx 0.05434$

# µ, sigma, monthly pr ratio, yearly pr ratio
vwrtn_basic_stats = apply(as.data.frame(lvwrtn_ts), 2, basicStats)
vwrtn_basic_stats$x['Mean',];
vwrtn_basic_stats$x['Stdev',];
mnthly_vw_pr_ratio = exp(vwrtn_basic_stats$x['Mean',]+vwrtn_basic_stats$x['Stdev',]^2/2);
mnthly_vw_pr_ratio; mnthly_vw_pr_ratio^12;

Sec. 2.5¶

See arima doc: link

MA(9) fitting with fixed coefs example (Sec. 2.5.3)¶

da = read.table("../AFTS_data/Ch02/m-ibm3dx2608.txt", header = T)
da[1:5,]

ewrtn_ts = ts(da$ewrtn, frequency = 12, start = c(1926, 1))
plot_time_fig(ts = ewrtn_ts, main = "Equal-weighted index return (Fig. 2.8)", xlab = "Year", ylab = "Simple Return")

par(bg = "white")
acf(da$ewrtn, main = "ACF: Equal-weighted index return (Fig. 2.8)")

find("arima")

Fit MA(9) process

ewrtn_ma_mod = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA), # Include the mean
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod; sqrt(ewrtn_ma_mod$sigma2)

Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA))

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22

ewrtn_ma_mod_1 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'ML',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_1; sqrt(ewrtn_ma_mod_1$sigma2)

Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "ML")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22

ewrtn_ma_mod_2 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'CSS',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_2; sqrt(ewrtn_ma_mod_2$sigma2)

Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "CSS")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1198    0    0    0    0    0  0.1232     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0313     0.0027

sigma^2 estimated as 0.005098:  part log likelihood = 1215.65

ewrtn_ma_mod_3 = arima(
    da$ewrtn, 
    order = c(0,0,9), 
    fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA),
    method = 'CSS-ML',
    transform.pars = F # Doesn't make difference
)
ewrtn_ma_mod_3; sqrt(ewrtn_ma_mod_3$sigma2)

Call:
arima(x = da$ewrtn, order = c(0, 0, 9), transform.pars = F, fixed = c(NA, 0, 
    NA, 0, 0, 0, 0, 0, NA, NA), method = "CSS-ML")

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1909    0  -0.1199    0    0    0    0    0  0.1227     0.0122
s.e.  0.0293    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005097:  log likelihood = 1215.61,  aic = -2421.22

Box-Ljung test (us 12 dof)

box_test_res = Box.test(ewrtn_ma_mod$residuals, lag = 12, type = 'Ljung');
box_test_res; 1-pchisq(box_test_res$statistic, df = 9); 1-pchisq(box_test_res$statistic, df = 12)

	Box-Ljung test

data:  ewrtn_ma_mod$residuals
X-squared = 17.604, df = 12, p-value = 0.1283

box_test_res_1 = Box.test(ewrtn_ma_mod_1$residuals, lag = 12, type = 'Ljung');
box_test_res_1; 1-pchisq(box_test_res_1$statistic, df = 9); 1-pchisq(box_test_res_1$statistic, df = 12)

	Box-Ljung test

data:  ewrtn_ma_mod_1$residuals
X-squared = 17.603, df = 12, p-value = 0.1283

box_test_res_2 = Box.test(ewrtn_ma_mod_2$residuals, lag = 12, type = 'Ljung');
box_test_res_2; 1-pchisq(box_test_res_2$statistic, df = 9); 1-pchisq(box_test_res_2$statistic, df = 12)

	Box-Ljung test

data:  ewrtn_ma_mod_2$residuals
X-squared = 17.621, df = 12, p-value = 0.1277

box_test_res_3 = Box.test(ewrtn_ma_mod_3$residuals, lag = 12, type = 'Ljung');
box_test_res_3; 1-pchisq(box_test_res_3$statistic, df = 9); 1-pchisq(box_test_res_3$statistic, df = 12)

	Box-Ljung test

data:  ewrtn_ma_mod_3$residuals
X-squared = 17.604, df = 12, p-value = 0.1283

MA(9) forecasting (Sec. 2.5.4)¶

Plot the forecasting

eotr = 986
h = 10
npts = 12
freq = 12
order = c(0,0,9)
fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA)
seasonal = NULL
ewrtn_fc_res = plot_arima_forecast_fig(
    da_ts=ewrtn_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=F,
    main="Forecasts from ARIMA(0,0,9) with non-zero mean for ewrtn (Fig. 2.7)", 
    xlab="Year", ylab="Return", ylim=c(-0.25, 0.25)
)

2007.083 ; 2009

Forecasting results

ewrtn_multistep_ahead_forecast_tb = comb_forecast_res(
    forecast_obj=ewrtn_fc_res, 
    da_ts=ewrtn_ts, 
    eotr=eotr, 
    freq=freq
)
ewrtn_multistep_ahead_forecast_tb

Forecast method: ARIMA(0,0,9) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, transform.pars = transform_pars, fixed = fixed, 
    method = method)

Coefficients:
         ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept
      0.1844    0  -0.1206    0    0    0    0    0  0.1218     0.0128
s.e.  0.0295    0   0.0338    0    0    0    0    0  0.0312     0.0027

sigma^2 estimated as 0.005066:  log likelihood = 1206.44,  aic = -2402.88

Error measures:
                        ME       RMSE        MAE      MPE     MAPE      MASE
Training set -1.916748e-05 0.07117457 0.04750136 94.10055 214.5367 0.6922456
                   ACF1
Training set 0.02241295

Forecasts:
         Point Forecast       Lo 80      Hi 80      Lo 95     Hi 95
Mar 2008    0.004284404 -0.08692947 0.09549828 -0.1352152 0.1437840
Apr 2008    0.013560426 -0.07919123 0.10631208 -0.1282910 0.1554118
May 2008    0.015025565 -0.07772609 0.10777722 -0.1268258 0.1568770
Jun 2008    0.014453940 -0.07894734 0.10785522 -0.1283910 0.1572989
Jul 2008    0.012047087 -0.08135419 0.10544836 -0.1307978 0.1548920
Aug 2008    0.001806804 -0.09159447 0.09520808 -0.1410381 0.1446517
Sep 2008    0.012211891 -0.08118938 0.10561317 -0.1306330 0.1550568
Oct 2008    0.005515984 -0.08788529 0.09891726 -0.1373289 0.1483609
Nov 2008    0.008514008 -0.08488727 0.10191528 -0.1343309 0.1513589
Dec 2008    0.012792573 -0.08126721 0.10685236 -0.1310595 0.1566446

Find Equal-weighted monthly and yearly expected return¶

# µ, sigma, monthly pr ratio, yearly pr ratio
lewrtn_ts = ts(log(1+ewrtn_ts), frequency = 12, start = c(1926, 1))
ewrtn_basic_stats = apply(as.data.frame(lewrtn_ts), 2, basicStats)
ewrtn_basic_stats$x['Mean',];
ewrtn_basic_stats$x['Stdev',];
mnthly_ew_pr_ratio = exp(ewrtn_basic_stats$x['Mean',]+ewrtn_basic_stats$x['Stdev',]^2/2);
mnthly_ew_pr_ratio; mnthly_ew_pr_ratio^12;

Sec. 2.6¶

3M log price (Sec. 2.6.3)¶

da = read.table("../AFTS_data/Ch02/m-3m4608.txt", header = T)
da[1:5,]

lrtn = log(1+da$rtn)
lrtn_ts = ts(lrtn, frequency = 12, start = c(1946, 2))

plot_time_fig(lrtn_ts, main = "3M stock (Fig. 2.9)", xlab = "Year", ylab = "log(Return)")

par(bg = "white")
acf(lrtn_ts, main = "3M stock (Fig. 2.9)")

Print EACF results (Table 2.5)
- eacf: link

find("eacf")

eacf_res <- eacf(lrtn_ts, ar.max = 6, ma.max = 12)
eacf_stats_tb = format(as.data.frame(eacf_res$eacf), digits = 3)
names(eacf_stats_tb) <- seq(from = 0, to = 12)
display(eacf_stats_tb)
display(eacf_res$symbol)
display(2/sqrt(length(lrtn_ts)))

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12
0 o o x o o x o o o x o  x  o 
1 x o x o o x o o o o o  x  o 
2 x x x o o x o o o o o  o  o 
3 x x x o o o o o o o o  o  o 
4 x o x o o o o o o o o  o  o 
5 x x x o x o o o o o o  o  o 
6 x x x x x o o o o o o  o  o

Find 3M monthly and yearly expected return¶

For the log(Return)
- $\mu_{monthly}\approx 0.010299$
- $\sigma_{monthly}\approx 0.063719$

basic_3m_stats <- apply(as.data.frame(lrtn), 2, basicStats)
basic_3m_stats
# µ, std, and std of mean
basic_3m_stats$lrtn['Mean',]; basic_3m_stats$lrtn['Stdev',]; basic_3m_stats$lrtn['Stdev',]/sqrt(basic_3m_stats$lrtn['nobs',])
# 1 + Monthly simple grow rate
mnthly_3m_pr_ratio = exp(basic_3m_stats$lrtn['Mean',]+basic_3m_stats$lrtn['Stdev',]^2/2) 
mnthly_3m_pr_ratio; mnthly_3m_pr_ratio^12

Sec. 2.7¶

3M log price¶

da = read.table("../AFTS_data/Ch02/m-3m4608.txt", header = T)
da[1:5,]

See the t-statistics and p-value of mean of log returns

lrtn = log(1+da$rtn)
lrtn_ts = ts(lrtn, frequency = 12, start = c(1946, 2))

t_stats = mean(lrtn)/(sd(lrtn)/sqrt(length(lrtn)))
pv = 2*(1-pnorm(abs(t_stats)))

mean(lrtn); 
sd(lrtn); 
sd(lrtn)/sqrt(length(lrtn)); # variance of the mean of the log returns
t_stats; 
pv

Construct two price series
- Using log return
- Using mean-corrected log return

lgPr1 = cumsum(lrtn)
lgPr2 = cumsum(lrtn-mean(lrtn))
lgPr1_ts = ts(lgPr1, frequency = 12, start = c(1946, 1))
lgPr2_ts = ts(lgPr2, frequency = 12, start = c(1946, 1))
ln_ts = ts(mean(lrtn)*seq(from = 1, to = length(lrtn), by = 1), frequency = 12, start = c(1946, 1))
par(bg = "white")
plot(ln_ts, col = 'black', type = 'l', main = '3M constructed log-pr series (Fig. 2.10)', xlab = "Year", ylab = "log(Price)")
# plot(lgPr1_ts, type = 'p', pch = 1, main = '3M constructed log-pr series', xlab = "Year", ylab = "Log Price")
points(lgPr1_ts, col = 'red', pch = 1)
points(lgPr2_ts, col = 'blue', pch = 3)
# Add legend
legend(
    "topleft", 
    legend = c("y=0.0103*t", "sum(r_t)", "sum(r_t-mean(t_t))"), 
    col = c("black", "red", "blue"), 
    lty = c(1, NA, NA),
    pch = c(NA, 1, 3)
)

See some basic stats for 3M return
- If we model the price series as a Geometric Brownian Motion, then the return series, i.e. diff(log(Price)), (with the constant time-step) should be sum of a constant mean and an iid Gaussian noise.
- The Geometric Brownian Motion can be seen as an unit-root process in time-series.

pr_ts = ts(exp(lgPr1), frequency = 12, start = c(1946, 1))
plot_time_fig(ts = pr_ts, main = "Price", xlab = "Year", ylab = "Price/Price_0")

Unit-root test (U.S. quarterly GDP)¶

da = read.table("../AFTS_data/Ch02/q-gdp4708.txt", header = T)
da[1:5,]

lgdp = log(da$gdp)
m1 = ar(diff(lgdp), method = 'mle')
m1$order

adfTest(lgdp, lags = 10, type = c("c"))

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 10
  STATISTIC:
    Dickey-Fuller: -1.6109
  P VALUE:
    0.4569 

Description:
 Sat Feb 17 20:09:22 2024 by user:

apply(as.data.frame(lgdp), 2, basicStats)

lgdp_ts = ts(lgdp, frequency = 4, start = c(1947 ,1))
plot_unit_root_figs(
    lgdp_ts, freq = 4, start = c(1947, 2), lag_max = 20,
    main = 'U.S. quarterly GDP (Fig. 2.11)', xlab = "Year", ylab = "ln(GDP)"
)

Full ARIMA workflow (S&P 500)¶

Unit-root testing
Order determination and parameter estimation
Box-Ljung test
Plot forcasting

# install.packages("RQuantLib")

also installing the dependencies ‘Rcpp’, ‘zoo’


Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done

load_quantlib_calendars(ql_calendars = NULL, from = 1950, to = 2008, financial = TRUE)

Calendar QuantLib/TARGET loaded

Calendar QuantLib/Argentina loaded

Calendar QuantLib/Australia loaded

Calendar QuantLib/Brazil loaded

Calendar QuantLib/Canada loaded

Calendar QuantLib/Canada/Settlement loaded

Calendar QuantLib/Canada/TSX loaded

Calendar QuantLib/China loaded

Calendar QuantLib/CzechRepublic loaded

Calendar QuantLib/Denmark loaded

Calendar QuantLib/Finland loaded

Calendar QuantLib/Germany loaded

Calendar QuantLib/Germany/FrankfurtStockExchange loaded

Calendar QuantLib/Germany/Settlement loaded

Calendar QuantLib/Germany/Xetra loaded

Calendar QuantLib/Germany/Eurex loaded

Calendar QuantLib/HongKong loaded

Calendar QuantLib/Hungary loaded

Calendar QuantLib/Iceland loaded

Calendar QuantLib/India loaded

Calendar QuantLib/Indonesia loaded

Calendar QuantLib/Italy loaded

Calendar QuantLib/Italy/Settlement loaded

Calendar QuantLib/Italy/Exchange loaded

Calendar QuantLib/Japan loaded

Calendar QuantLib/Mexico loaded

Calendar QuantLib/NewZealand loaded

Calendar QuantLib/Norway loaded

Calendar QuantLib/Poland loaded

Calendar QuantLib/Russia loaded

Calendar QuantLib/SaudiArabia loaded

Calendar QuantLib/Singapore loaded

Calendar QuantLib/Slovakia loaded

Calendar QuantLib/SouthAfrica loaded

Calendar QuantLib/SouthKorea loaded

Calendar QuantLib/SouthKorea/KRX loaded

Calendar QuantLib/Sweden loaded

Calendar QuantLib/Switzerland loaded

Calendar QuantLib/Taiwan loaded

Calendar QuantLib/Turkey loaded

Calendar QuantLib/Ukraine loaded

Calendar QuantLib/UnitedKingdom loaded

Calendar QuantLib/UnitedKingdom/Settlement loaded

Calendar QuantLib/UnitedKingdom/Exchange loaded

Calendar QuantLib/UnitedKingdom/Metals loaded

Calendar QuantLib/UnitedStates loaded

Calendar QuantLib/UnitedStates/Settlement loaded

Calendar QuantLib/UnitedStates/NYSE loaded

Calendar QuantLib/UnitedStates/GovernmentBond loaded

Calendar QuantLib/UnitedStates/NERC loaded

da = read.table("../AFTS_data/Ch02/d-sp55008.txt", header = T)
da[1:10,]; dim(da)

Unit-root testing¶

lsp5 = log(da$close)
par(bg = "white")
lsp5_ts = ts(lsp5, frequency = 252, start = c(1950 , 1, 3))
plot(lsp5_ts, col = 'black', type = 'l', main = 'S&P 500 (Fig. 2.12)', xlab = "Year", ylab = "log(S&P 500)")

m2 = ar(diff(lsp5), method = 'mle')
m2$order

adfTest(lsp5, lags = 2, type = ("ct"))

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 2
  STATISTIC:
    Dickey-Fuller: -2.0179
  P VALUE:
    0.5708 

Description:
 Sat Feb 17 20:49:04 2024 by user:

log(length(lsp5))

adfTest(lsp5, lags = 15, type = ("ct"))

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 15
  STATISTIC:
    Dickey-Fuller: -1.9946
  P VALUE:
    0.5807 

Description:
 Wed Feb 14 03:52:33 2024 by user:

Order determination and parameter estimations¶

- With one-differencing

diff_lsp5_ts = ts(diff(lsp5), frequency = 252, start = c(1950 , 1, 4))
length(lsp5); length(diff_lsp5_ts);
display(head(lsp5, 10))
display(head(diff_lsp5_ts, 10))

plot_unit_root_figs(
    lsp5_ts, freq = 252, start = c(1950 , 1, 4), lag_max = 20,
    main = 'S&P500 Index', xlab = "Year", ylab = "ln(S&P500)"
)

diff_lsp5_ts = ts(diff(lsp5), frequency = 252, start = c(1950 , 1, 4))

par(bg = 'white')
pacf(diff_lsp5_ts)

par(bg = 'white')
acf(diff_lsp5_ts)

sp5_eacf_res <- eacf(diff_lsp5_ts, ar.max = 6, ma.max = 12)
sp5_eacf_stats_tb = format(as.data.frame(sp5_eacf_res$eacf), digits = 3)
names(sp5_eacf_stats_tb) <- seq(from = 0, to = 12)
display(sp5_eacf_stats_tb)
display(sp5_eacf_res$symbol)
# pp67, asymptotic standard error of EACF
display(2/sqrt(length(diff_lsp5_ts)))

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12
0 x x o o o x o o o o o  x  o 
1 x x o o o o o o o o o  x  o 
2 x o o o o o o o o o o  x  o 
3 x x o o o o o o o o o  x  o 
4 x o o x o o o o o o o  x  o 
5 o x x x x x o o o o o  o  o 
6 x x x x x x o o o o o  o  o

diff_sp5_mod = arima(diff_lsp5_ts, order = c(1,0,2), include.mean = F)
diff_sp5_mod

Call:
arima(x = diff_lsp5_ts, order = c(1, 0, 2), include.mean = F)

Coefficients:
         ar1     ma1      ma2
      0.0726  0.0001  -0.0373
s.e.  0.2893  0.2894   0.0251

sigma^2 estimated as 8.077e-05:  log likelihood = 48279.03,  aic = -96552.06

sp5_mod = arima(lsp5_ts, order = c(1,1,2), include.mean = T)
sp5_mod

Call:
arima(x = lsp5_ts, order = c(1, 1, 2), include.mean = T)

Coefficients:
         ar1     ma1      ma2
      0.0726  0.0001  -0.0373
s.e.  0.2893  0.2894   0.0251

sigma^2 estimated as 8.077e-05:  log likelihood = 48279.03,  aic = -96552.06

Box-Ljung test¶

pp33 very briefly talks about how to choose lag in the Box.test
Null hypothesis cannot be reject, so that the model is adequte.

diff_sp5_box_test = Box.test(diff_sp5_mod$residuals, lag = 5, type = 'Ljung')
diff_sp5_box_test

	Box-Ljung test

data:  diff_sp5_mod$residuals
X-squared = 1.6681, df = 5, p-value = 0.8929

diff_sp5_box_test_1 = Box.test(diff_sp5_mod$residuals, lag = 12, type = 'Ljung')
diff_sp5_box_test_1

	Box-Ljung test

data:  diff_sp5_mod$residuals
X-squared = 20.642, df = 12, p-value = 0.05587

sp5_box_test = Box.test(sp5_mod$residuals, lag = 5, type = 'Ljung')
sp5_box_test

	Box-Ljung test

data:  sp5_mod$residuals
X-squared = 1.6646, df = 5, p-value = 0.8933

sp5_box_test_1 = Box.test(sp5_mod$residuals, lag = 12, type = 'Ljung')
sp5_box_test_1

	Box-Ljung test

data:  sp5_mod$residuals
X-squared = 20.656, df = 12, p-value = 0.05565

Plot forcasting¶

diff_sp5_mod$nobs

npts = 22
eotr = 14661-npts
h = npts
freq = 252
order_1 = c(1,0,2)
fixed = NULL
seasonal = NULL
diff_sp5_fc_res = plot_arima_forecast_fig(
    da_ts=diff_lsp5_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order_1, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=TRUE,
    main="Forecasts from ARIMA(1,0,2) with\n non-zero mean for ln(Return) of s&p500", 
    xlab="Year", ylab="ln(Return)", ylim=c(-0.05, 0.05)
)

[1] 14639
2008 ; 2008.179

order_2 = c(1,1,2)
sp5_fc_res = plot_arima_forecast_fig(
    da_ts=lsp5_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order_2, seasonal=seasonal, fixed=fixed, method='ML', transform.pars=TRUE,
    main="Forecasts from ARIMA(1,1,2) with\n non-zero mean for ln(Price) of s&p500", 
    xlab="Year", ylab="ln(Return)", ylim=c(7, 7.3)
)

[1] 14638
2008 ; 2008.183

Find S&P500 monthly and yearly expected return¶

We can think of
For the log(Return)
- $\mu_{daily}\approx 0.000299$
- $\sigma_{daily}\approx 0.009011$

# µ, std, and std of mean; daily, monthly, yearly pr-ratio
sp5_basic_stats <- apply(as.data.frame(diff(lsp5)), 2, basicStats)
display(sp5_basic_stats)
daily_exp_pr_ratio = exp(sp5_basic_stats$`diff(lsp5)`['Mean',]+sp5_basic_stats$`diff(lsp5)`['Stdev',]^2/2)
mnthly_exp_pr_ratio = daily_exp_pr_ratio^21
yrly_exp_pr_ratio = daily_exp_pr_ratio^252
sp5_basic_stats$`diff(lsp5)`['Mean',]; 
sp5_basic_stats$`diff(lsp5)`['Stdev',];
daily_exp_pr_ratio; mnthly_exp_pr_ratio; yrly_exp_pr_ratio;

Sec. 2.8¶

J&J quarterly earnings¶

Starting at 1960 first quarter.

da = read.table("../AFTS_data/Ch02/q-jnj.txt", header = F) # needs to set `header = F`
head(da); length(da$V1)

Plot the acf of JnJ quarterly earnings¶

freq = 4
jnj_er = da$V1
jnj_er_ts = ts(jnj_er, frequency = freq, start = c(1960, 1))
jnj_lg_er = log(jnj_er)
jnj_lg_er_ts = ts(jnj_lg_er, frequency = freq, start = c(1960, 1))
jnj_lg_er_ts

plot_time_fig(jnj_er_ts, main = "JnJ quarterly earnings (Fig. 2.13 (a))", xlab = "Year", ylab = "Earning")
plot_time_fig(jnj_lg_er_ts, main = "JnJ quarterly earnings (Fig. 2.13 (b))", xlab = "Year", ylab = "log(Earning)")

plot_pacf_acf(jnj_lg_er, main = "lg(earning) (Fig. 2.14 (a))")
plot_pacf_acf(diff(jnj_lg_er), main = "diff(lg(earning)) (Fig. 2.14 (b))")
plot_pacf_acf(diff(jnj_lg_er, lag = freq), main = "diff(lg(earning), 4) (Fig. 2.14 (c))")
plot_pacf_acf(diff(diff(jnj_lg_er), lag = freq), main = "diff(diff(lg(earning)), 4) (Fig. 2.14 (d))")
plot_pacf_acf(diff(diff(jnj_lg_er, lag = freq)), main = "diff(diff(lg(earning), 4)) (Fig. 2.14 (d)\')")

Apply airline model to the JnJ data.¶

find('arima')

multi_seas_mod_1 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'ML',
    include.mean = F
)
summary(multi_seas_mod_1); sqrt(multi_seas_mod_1$sigma2)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "ML")

Coefficients:
          ma1     sma1
      -0.6809  -0.3146
s.e.   0.0982   0.1070

sigma^2 estimated as 0.007931:  log likelihood = 78.38,  aic = -152.75

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

multi_seas_mod_2 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'CSS-ML',
    include.mean = F
)
summary(multi_seas_mod_2); sqrt(multi_seas_mod_2$sigma2)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "CSS-ML")

Coefficients:
          ma1     sma1
      -0.6809  -0.3146
s.e.   0.0982   0.1070

sigma^2 estimated as 0.007931:  log likelihood = 78.38,  aic = -152.75

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

multi_seas_mod_3 <- arima(
    x = diff(diff(jnj_lg_er_ts), lag = freq),
    order = c(0,0,1),
    seasonal = list(order = c(0,0,1), period = freq),
    method = 'CSS',
    include.mean = F
)
summary(multi_seas_mod_3); sqrt(multi_seas_mod_3$sigma2)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(0, 0, 1), seasonal = list(order = c(0, 
    0, 1), period = freq), include.mean = F, method = "CSS")

Coefficients:
          ma1     sma1
      -0.5228  -0.2752
s.e.   0.0956   0.0992

sigma^2 estimated as 0.008665:  part log likelihood = 75.47

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

# multi_seas_mod_4 <- arima(
#     x = diff(diff(jnj_lg_er_ts), lag = freq),
#     order = c(1,0,1),
#     seasonal = list(order = c(1,0,1), period = freq),
#     method = 'ML',
#     include.mean = F,
#     transform.pars = F
# )
# # multi_seas_mod_4 <- arima(
# #     x = diff(diff(jnj_lg_er_ts), lag = freq),
# #     order = c(1,0,1),
# #     seasonal = list(order = c(1,0,1), period = freq),
# #     fixed = c(1, NA, 1, NA), 
# #     method = 'ML',
# #     include.mean = F,
# #     transform.pars = F
# # )
# summary(multi_seas_mod_4); sqrt(multi_seas_mod_4$sigma2)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = diff(diff(jnj_lg_er_ts), lag = freq), order = c(1, 0, 1), seasonal = list(order = c(1, 
    0, 1), period = freq), include.mean = F, transform.pars = F, method = "ML")

Coefficients:
          ar1      ma1     sar1     sma1
      -0.0180  -0.6626  -0.2075  -0.1382
s.e.   0.2216   0.1815   0.3080   0.2950

sigma^2 estimated as 0.007889:  log likelihood = 78.57,  aic = -149.14

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

Box-Ljung test¶

Cannot reject the null-hypothesis. So the model appears to be adequate.

box_test <- Box.test(multi_seas_mod_1$residuals, lag = 12, type = 'Ljung')
box_test

	Box-Ljung test

data:  multi_seas_mod_1$residuals
X-squared = 10.15, df = 12, p-value = 0.6028

Forecast¶

(0):

$$ \begin{eqnarray} z_t &=& (1-B)(1-B^4)x_t = (1-\theta B)(1-\beta B^4)a_t = a_t-\theta a_{t-1}-\beta a_{t-4}+\theta \beta a_{t-5} \\ y_t &=& (1-B^4)x_t \\ z_t &=& (1-B)y_t \end{eqnarray} $$

$$ \begin{eqnarray} y_{t+1}&=&z_{t+1}+y_t \\ x_{t+4}&=&y_{t+4}+x_t \end{eqnarray} $$

or

$$ \begin{eqnarray} y_{t}&=&z_{t}+y_{t-1} \\ x_{t}&=&y_{t}+x_{t-4} \end{eqnarray} $$

The expectation is

$$ \begin{eqnarray} E[y_{t+1}]&=&E[z_{t+1}]+E[y_t] \\ E[x_{t+4}]&=&E[y_{t+4}]+E[x_t] \end{eqnarray} $$

(1) (following are stationary results, but for forecasting needs to consider the time-zero):

$$ \begin{eqnarray} Var[z_t]&=&(1+\theta^2+\beta^2+\theta^2\beta^2)\sigma^2 \\ Cov(z_{t+1},z_t)&=&-\theta(1+\beta^2)\sigma^2 \\ Cov(z_{t+2},z_t)&=&0 \\ Cov(z_{t+3},z_t)&=&\theta\beta\sigma^2 \\ Cov(z_{t+4},z_t)&=&-\beta(1+\theta^2)\sigma^2 \\ Cov(z_{t+5},z_t)&=&\theta\beta\sigma^2 \\ Cov(z_{t+i},z_t)&=&0, \forall i>5 \end{eqnarray} $$

(2):

$$ \begin{eqnarray} Var[y_t]&=&Var[z_t]+Var[y_{t-1}]+Cov(z_t, y_{t-1}) \\ Cov(y_{t+1}, y_t)&=&Cov(z_{t+1}, y_t)+Var[y_t]=Cov(z_{t+1},z_t+z_{t-1}+z_{t-2}+z_{t-3}+z_{t-4})+Var[y_t] \\ Cov(y_{t+2}, y_t)&=&Cov(z_{t+2}, y_t)+Cov(y_{t+1}, y_t)=Cov(z_{t+2},z_t+z_{t-1}+z_{t-2}+z_{t-3})+Cov(y_{t+1}, y_t) \\ Cov(y_{t+3}, y_t)&=&Cov(z_{t+3}, y_t)+Cov(y_{t+2}, y_t)=Cov(z_{t+3},z_t+z_{t-1}+z_{t-2})+Cov(y_{t+2}, y_t) \\ Cov(y_{t+4}, y_t)&=&Cov(z_{t+4}, y_t)+Cov(y_{t+3}, y_t)=Cov(z_{t+4},z_t+z_{t-1})+Cov(y_{t+3}, y_t) \\ Cov(y_{t+5}, y_t)&=&Cov(z_{t+5}, y_t)+Cov(y_{t+4}, y_t)=Cov(z_{t+5},z_t)+Cov(y_{t+4}, y_t) \\ Cov(y_{t+i}, y_t)&=&Cov(y_{t+5}, y_t), \forall i>5 \end{eqnarray} $$

(3):

$$ \begin{eqnarray} Var[x_t]&=&Var[y_t]+Var[x_{t-4}]+Cov(y_t,y_{t-4}+y_{t-8}+\cdots) \end{eqnarray} $$

# plot_arima_forecast_fig <- function(
#     da_ts, eotr, h, npts, frequency,
#     order, seasonal, fixed, method, include.mean, transform.pars,
#     main, xlab, ylab, ylim = NULL) {
#     par(bg = "white")
#     # arima model
#     tr_da_ts <- ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
#     if (is.null(transform.pars)) {
#         ts_fm3 <- arima(tr_da_ts, order = order, fixed = fixed, seasonal = seasonal, method = method, include.mean = include.mean)
#     } else {
#         ts_fm3 <- arima(tr_da_ts, order = order, fixed = fixed, seasonal = seasonal, method = method, include.mean = include.mean, transform.pars = transform.pars)
#     }
#     print(ts_fm3$nobs)
#     # Forecast
#     ts_fm3$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
#     ts_fc_res <- forecast(ts_fm3, h = h)
#     # Plot forecast
#     if (is.null(npts)) {
#         npts <- eotr
#     }
#     xmin <- time(da_ts)[eotr] - npts / frequency
#     xmax <- time(da_ts)[eotr] + (max(h, length(da_ts) - eotr) + 1) / frequency
#     cat(xmin, ";", xmax)

#     plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
#     # Plot forecast mean (prepend the last observed data in the training dataset)
#     dummy_1st_fmean_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_fmean_ts)
#     points(dummy_1st_fmean_ts, pch = 1)
#     # Plot confidence interval (95%)
#     dummy_1st_flower_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[, 2])), frequency = frequency, start = end(tr_da_ts))
#     dummy_1st_fupper_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[, 2])), frequency = frequency, start = end(tr_da_ts))
#     lines(dummy_1st_flower_ts, lty = 2)
#     lines(dummy_1st_fupper_ts, lty = 2)
#     # Plot original data
#     orig_plot_ts <- ts(da_ts[(eotr - npts + 1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr] - (npts - 1) / frequency)
#     lines(orig_plot_ts)
#     points(orig_plot_ts, pch = 19)
#     ts_fc_res
# }

# comb_forecast_res <- function(forecast_obj, da_ts, eotr, freq) {
#     display(summary(forecast_obj))
#     fc_std_err = (forecast_obj$upper[,2]-forecast_obj$lower[,2])/2/qnorm(p = 0.975)
#     actual_ts = ts(da_ts[(eotr+1):length(da_ts)], frequency = freq, start = time(da_ts)[eotr+1])
#     display(forecast_obj$mean); display(fc_std_err); display(actual_ts)
#     multistep_ahead_forecast_tb = cbind(forecast_obj$mean, fc_std_err, actual_ts)
#     dimnames(multistep_ahead_forecast_tb)[[2]] <- c("Forecast", "Std. Error", "Actual")
#     multistep_ahead_forecast_tb
# }

npts = 6
h = 8
eotr = 79-h
freq = 4
order = c(0,0,1)
seasonal = list(order = order, period = freq)
fixed = NULL
include.mean = F
diff_seasonal_reg_jnj_ts = ts(diff(diff(jnj_lg_er_ts), lag = freq), frequency = freq, start = c(1961, 2))
diff_seasonal_reg_jnj_fc_res = plot_arima_forecast_fig(
    da_ts=diff_seasonal_reg_jnj_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', include.mean=include.mean, transform.pars=TRUE,
    main="Forecasts from ARIMA(0,0,1)(0,0,1)[4] with\n non-zero mean for (1-B)(1-B^4)ln(QtlyEarning) of JnJ", 
    xlab="Year", ylab="(1-B)(1-B^4)ln(QtlyEarning)"#, ylim=c(-0.05, 0.05)
)
summary(diff_seasonal_reg_jnj_fc_res)

[1] 71
1977.25 ; 1981

Forecast method: ARIMA(0,0,1)(0,0,1)[4] with zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME      RMSE        MAE  MPE MAPE      MASE       ACF1
Training set 0.01083226 0.0916983 0.07697376 -Inf  Inf 0.5682199 0.05218151

Forecasts:
        Point Forecast       Lo 80     Hi 80       Lo 95     Hi 95
1979 Q1     0.09037218 -0.02714392 0.2078883 -0.08935319 0.2700975
1979 Q2     0.02639844 -0.11413996 0.1669369 -0.18853650 0.2413334
1979 Q3    -0.02343313 -0.16397154 0.1171053 -0.23836807 0.1915018
1979 Q4     0.06359381 -0.07694460 0.2041322 -0.15134113 0.2785288
1980 Q1    -0.03378421 -0.18019118 0.1126228 -0.25769434 0.1901259
1980 Q2     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q3     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q4     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624

jnj_multistep_ahead_forecast_tb <- comb_forecast_res(diff_seasonal_reg_jnj_fc_res, diff_seasonal_reg_jnj_ts, eotr, freq)
jnj_multistep_ahead_forecast_tb

Forecast method: ARIMA(0,0,1)(0,0,1)[4] with zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME      RMSE        MAE  MPE MAPE      MASE       ACF1
Training set 0.01083226 0.0916983 0.07697376 -Inf  Inf 0.5682199 0.05218151

Forecasts:
        Point Forecast       Lo 80     Hi 80       Lo 95     Hi 95
1979 Q1     0.09037218 -0.02714392 0.2078883 -0.08935319 0.2700975
1979 Q2     0.02639844 -0.11413996 0.1669369 -0.18853650 0.2413334
1979 Q3    -0.02343313 -0.16397154 0.1171053 -0.23836807 0.1915018
1979 Q4     0.06359381 -0.07694460 0.2041322 -0.15134113 0.2785288
1980 Q1    -0.03378421 -0.18019118 0.1126228 -0.25769434 0.1901259
1980 Q2     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q3     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624
1980 Q4     0.00000000 -0.14886044 0.1488604 -0.22766240 0.2276624

Use DP to calculate forecasting expectation and variance. $$ \begin{eqnarray} E[y_{t+1}]&=&E[z_{t+1}]+E[y_t] \\ E[x_{t+4}]&=&E[y_{t+4}]+E[x_t] \end{eqnarray} $$

# Assume freq=s>=4
# z_t
diff_seasonal_reg_jnj_ts = ts(diff(diff(jnj_lg_er_ts), lag = freq), frequency = freq, start = c(1961, 2))
# y_t
diff_seasonal_jnj_ts = ts(diff(jnj_lg_er_ts, lag = freq), frequency = freq, start = c(1961, 1))
fstart = c(1979, 1)
iz = 71
iy = 72
ix = 76
h = 8
# E[y_t]
e_y <- list()
for (i in 1:h) {
    if (i==1) {
        e_y[[paste(i)]] = diff_seasonal_reg_jnj_fc_res$mean[i] + diff_seasonal_jnj_ts[length(diff_seasonal_jnj_ts)-h]
    } else {
        e_y[[paste(i)]] = diff_seasonal_reg_jnj_fc_res$mean[i] + e_y[[paste(i-1)]]
    }
}
# E[x_t]
e_x <- list()
for (i in 1:h) {
    if (i<=freq) {
        e_x[[paste(i)]] = e_y[[paste(i)]] + jnj_lg_er_ts[length(jnj_lg_er_ts)-h+i-freq]
    } else {
        e_x[[paste(i)]] = e_y[[paste(i)]] + e_x[[paste(i-freq)]]
    }
}
theta = -as.numeric(diff_seasonal_reg_jnj_fc_res$model$coef[1])
beta = -as.numeric(diff_seasonal_reg_jnj_fc_res$model$coef[2])
sigma2 = as.numeric(diff_seasonal_reg_jnj_fc_res$model$sigma2)
gen_key <- function(i,j) {paste("(",i,",",j,")", sep = "")}
# V[z_t]
# The following are unconditional, we need to calculate conditional
# v_z <- list(
#     "0"=as.numeric((1+theta^2+beta^2+theta^2*beta^2)*sigma2), 
#     "1"=as.numeric(-theta*(1+beta^2)*sigma2),
#     "2"=0,
#     "3"=as.numeric(theta*beta*sigma2),
#     "4"=as.numeric(-beta*(1+theta^2)*sigma2),
#     "5"=as.numeric(theta*beta*sigma2)
# )
v_z <- list()
for (i in 1:h) {
    for (j in 0:h) {
        v_z[[gen_key(i,j)]] = 0
        if (j==0) {
            v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + sigma2
            if (i-1>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + theta^2*sigma2}
            if (i-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + beta^2*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] + theta^2*beta^2*sigma2}
        } else if (j==1) {
            v_z[[gen_key(i,j)]] = -theta*sigma2
            if (i-freq>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] - theta*beta^2*sigma2}
        } else if (j==freq-1) {
            if (i-1>0) {v_z[[gen_key(i,j)]] = theta*beta*sigma2}
        } else if (j==freq) {
            v_z[[gen_key(i,j)]] = -beta*sigma2
            if (i-1>0) {v_z[[gen_key(i,j)]] = v_z[[gen_key(i,j)]] - beta*theta^2*sigma2}
        } else if (j==freq+1) {v_z[[gen_key(i,j)]] = theta*beta*sigma2}
        else {v_z[[gen_key(i,j)]] = 0}
    }
    for (j in 1:h) {
        v_z[[gen_key(i,-j)]] = 0
        if (j==1) {
            if (i-1>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - theta*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - theta*beta^2*sigma2}
        } else if (j==freq-1) {
            if (i-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] + theta*beta*sigma2}
        } else if (j==freq) {
            if (i-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - beta*sigma2}
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] - beta*theta^2*sigma2}
        } else if (j==freq+1) {
            if (i-1-freq>0) {v_z[[gen_key(i,-j)]] = v_z[[gen_key(i,-j)]] + theta*beta*sigma2}
        }
    }
}
# V[y_t]
v_tau <- list(
    "0"=sigma2,
    "1"=-theta*sigma2,
    "4"=-beta*sigma2,
    "5"=theta*beta*sigma2
)
v_y <- list()
for (i in 1:h) {
    v_y[[gen_key(i,0)]] = v_z[[gen_key(i,0)]]
    if (i-1>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] + v_y[[gen_key(i-1,0)]] - theta*v_tau[["0"]]}
    if (i-freq>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] - beta*(v_tau[["0"]]+v_tau[["1"]])}
    if (i-1-freq>0) {v_y[[gen_key(i,0)]] = v_y[[gen_key(i,0)]] + theta*beta*(v_tau[["0"]]+v_tau[["1"]]+v_tau[["4"]])}
    for (j in 1:h) {
        v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j-1)]]
        if (j<=freq+1) {
            for (k in j:(freq+1)) {
                if (i+j-k>0) {v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j)]] + v_z[[gen_key(i+j-k,k)]]}
                else {v_y[[gen_key(i,j)]] = v_y[[gen_key(i,j)]] + v_z[[gen_key(i+j,-k)]]}
            }
        }
    }
}
# V[x_t]
v_x <- list()
for (i in 1:h) {
    v_x[[paste(i)]] = v_y[[gen_key(i,0)]]
    if (i>freq) {v_x[[paste(i)]] = v_x[[paste(i)]] + v_x[[paste(i-freq)]]}
    j = freq
    while (TRUE) {
        if (i>j) {
            v_x[[paste(i)]] = v_x[[paste(i)]] + v_y[[gen_key(i-j,j)]]
            j = j + freq
        }
        else {break}
    }
}

e_y; e_x

# My calculation of std.error match the results from the package
for (i in 1:h) {
    cat(i, ":", sqrt(v_z[[gen_key(i,0)]]), "\n")
}
jnj_multistep_ahead_forecast_tb

1 : 0.0916983 
2 : 0.1096627 
3 : 0.1096627 
4 : 0.1096627 
5 : 0.114242 
6 : 0.1161564 
7 : 0.1161564 
8 : 0.1161564

v_z

v_y

v_x

theta; beta; sigma2

Plot the forecasting
- The resulting plot is somewhat different from the textbook Fig. 2.15.

e_qer_no_var_arr <- c()
e_lg_qer_arr <- c()
e_qer_arr <- c()
v_qer_arr <- c()
for (i in 1:h) {
    e_qer_no_var_arr <- c(e_qer_no_var_arr, exp(e_x[[paste(i)]]))
    e_lg_qer_arr <- c(e_lg_qer_arr, e_x[[paste(i)]])
    e_qer_arr <- c(e_qer_arr, exp(e_x[[paste(i)]]+v_x[[paste(i)]]/2))
    v_qer_arr <- c(v_qer_arr, v_x[[paste(i)]])
    cat(i, ":", 
        exp(e_x[[paste(i)]]), exp(e_x[[paste(i)]]+v_x[[paste(i)]]/2), 
        e_x[[paste(i)]], v_x[[paste(i)]], "\n",
        sep = ","
    )
}

1,:,13.27175,13.32766,2.585637,0.008408578,
2,:,13.83323,13.93681,2.627074,0.01491936,
3,:,13.61369,13.76034,2.611076,0.02143014,
4,:,10.63887,10.78855,2.364515,0.02794093,
5,:,15.32053,15.65217,2.729194,0.04283245,
6,:,15.96869,16.47529,2.77063,0.06246335,
7,:,15.71525,16.37374,2.754632,0.08209424,
8,:,12.28121,12.92202,2.508071,0.1017251,

help(legend)

par(bg = "white")
da_ts = jnj_er_ts
lg_da_ts = jnj_lg_er_ts
eotr = length(da_ts)-h
tr_da_ts <- ts(da_ts[1:eotr], frequency = freq, start = start(da_ts))
xmin <- time(da_ts)[eotr] - (npts+1) / freq
xmax <- time(da_ts)[eotr] + (h+1) / freq
xlim = c(xmin, xmax)
ylim = c(0, 30)
cat(xmin, ";", xmax)
# Label 1: Actual Earning Line
plot(da_ts, 
     xlim = xlim,
     ylim = ylim,
     main = "JnJ QtlyEarning Multiplicative\nSeasonal Model Forecasting (Fig. 2.15)", 
     xlab = "Year", ylab = "QtlyEarning", lty=1
)
# Plot forecast mean (prepend the last observed data in the training dataset)
dummy_1st_fmean_ts <- ts(
    c(c(da_ts[eotr]), as.numeric(e_qer_arr)), 
    frequency = freq, start = end(tr_da_ts)
)
# Label 2: NULL
lines(dummy_1st_fmean_ts)
# Label 3: Forecast Mean
points(dummy_1st_fmean_ts, pch = 1)
# Plot confidence interval (95%)
t_stats = qnorm(p = 0.975, mean = 0, sd = 1)
dummy_1st_flower_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), e_lg_qer_arr)-t_stats*c(c(0), sqrt(v_qer_arr))), # lower CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
dummy_1st_fupper_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), e_lg_qer_arr)+t_stats*c(c(0), sqrt(v_qer_arr))), # upper CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
# Label 4: Forecast 95% Lower Bound 
lines(dummy_1st_flower_ts, lty = 2)
# Label 5: Forecast 95% Upper Bound
lines(dummy_1st_fupper_ts, lty = 2)
# Plot original data
orig_plot_ts <- ts(
    da_ts[(eotr - npts + 1):length(da_ts)], 
    frequency = freq, start = time(da_ts)[eotr] - (npts - 1) / freq
)
# Label 6: NULL
lines(orig_plot_ts)
# Label 7: Actual Earning
points(orig_plot_ts, pch = 19)
legend(
    "topleft", 
    legend = c(
        "Actual Earning Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
        "Forecast 95% Upper Bound", NULL, "Actual Earning"
    ), 
#     col = c("black", "red", "blue"), 
    lty = c(1, NA, 2, 2, NA),
    pch = c(NA, 1, NA, NA, 19)
)

1977 ; 1981

Forecast: Directly apply multiplicative seasonal model¶

The forecast mean from my calculation and the R package are close to each other.
This is gives much closer CI to the textbook results than my calculation. My calculation inflate the forecasting variance.

length(jnj_lg_er_ts)

npts = 6
h = 8
eotr = length(jnj_lg_er_ts)-h
freq = 4
order = c(0,1,1) # From the solution of 2-6
seasonal = list(order = order, period = freq)
fixed = NULL
include.mean = F
multi_seas_mod_jnj_ts = ts(jnj_lg_er_ts, frequency = freq, start = c(1960, 1))
multi_seas_mod_jnj_fc_res = plot_arima_forecast_fig(
    da_ts=multi_seas_mod_jnj_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', include.mean=include.mean, transform.pars=TRUE,
    main="Forecasts from ARIMA(0,1,1)(0,1,1)[4] with non-zero mean for\n(1-B)(1-B^4)ln(QtlyEarning) = (1-theta*B)(1-Theta*B^4)a_t of JnJ", 
    xlab="Year", ylab="QtlyEarning", ylim=c(2, 3.2)
)
summary(multi_seas_mod_jnj_fc_res)

[1] 71
1977.25 ; 1981

Forecast method: ARIMA(0,1,1)(0,1,1)[4]

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set 0.01012177 0.08863074 0.07193184 NaN  Inf 0.4398973 0.04998853

Forecasts:
        Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1979 Q1       2.585637 2.468121 2.703154 2.405912 2.765363
1979 Q2       2.627074 2.502794 2.751353 2.437005 2.817142
1979 Q3       2.611076 2.480383 2.741768 2.411198 2.810953
1979 Q4       2.364514 2.227709 2.501320 2.155288 2.573741
1980 Q1       2.729193 2.549233 2.909154 2.453968 3.004419
1980 Q2       2.770630 2.578687 2.962573 2.477078 3.064181
1980 Q3       2.754632 2.551411 2.957852 2.443833 3.065430
1980 Q4       2.508070 2.294167 2.721974 2.180933 2.835208

multi_seas_mod_jnj_fc_tb = comb_forecast_res(
    multi_seas_mod_jnj_fc_res, 
    da_ts = multi_seas_mod_jnj_ts,
    eotr=eotr,
    freq=freq
)
multi_seas_mod_jnj_fc_df = as.data.frame(multi_seas_mod_jnj_fc_tb)
multi_seas_mod_jnj_fc_tb

Forecast method: ARIMA(0,1,1)(0,1,1)[4]

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.6559  -0.3492
s.e.   0.1094   0.1104

sigma^2 estimated as 0.008409:  log likelihood = 68.28,  aic = -132.57

Error measures:
                     ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set 0.01012177 0.08863074 0.07193184 NaN  Inf 0.4398973 0.04998853

Forecasts:
        Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1979 Q1       2.585637 2.468121 2.703154 2.405912 2.765363
1979 Q2       2.627074 2.502794 2.751353 2.437005 2.817142
1979 Q3       2.611076 2.480383 2.741768 2.411198 2.810953
1979 Q4       2.364514 2.227709 2.501320 2.155288 2.573741
1980 Q1       2.729193 2.549233 2.909154 2.453968 3.004419
1980 Q2       2.770630 2.578687 2.962573 2.477078 3.064181
1980 Q3       2.754632 2.551411 2.957852 2.443833 3.065430
1980 Q4       2.508070 2.294167 2.721974 2.180933 2.835208

multi_seas_mod_jnj_fc_df$Forecast

As comparison, the following are expectation from my calculation (my calculation is close to the result from package)

$`1`
2.58563748050113
$`2`
2.62707380269453
$`3`
2.61107565060245
$`4`
2.36451453244944
$`5`
2.72919356620702
$`6`
2.77062988840042
$`7`
2.75463173630834
$`8`
2.50807061815533

multi_seas_mod_jnj_fc_df$`Std. Error`^2

As comparison, the following are variance from my calculation (my calculation somhow inflate the variance)

$`1`
0.00840857826352835
$`2`
0.0149193615519221
$`3`
0.0214301448403158
$`4`
0.0279409281287095
$`5`
0.0428324475174856
$`6`
0.0624633460377802
$`7`
0.0820942445580749
$`8`
0.101725143078369

par(bg = "white")
da_ts = jnj_er_ts
lg_da_ts = jnj_lg_er_ts
eotr = length(da_ts)-h
tr_da_ts <- ts(da_ts[1:eotr], frequency = freq, start = start(da_ts))
xmin <- time(da_ts)[eotr] - (npts+1) / freq
xmax <- time(da_ts)[eotr] + (h+1) / freq
xlim = c(xmin, xmax)
ylim = c(0, 30)
cat(xmin, ";", xmax)
# Label 1: Actual Earning Line
plot(da_ts, 
     xlim = xlim,
     ylim = ylim,
     main = "JnJ QtlyEarning Multiplicative\nSeasonal Model Forecasting (Fig. 2.15)", 
     xlab = "Year", ylab = "QtlyEarning", lty=1
)
# Plot forecast mean (prepend the last observed data in the training dataset)
multi_seas_e_qer_arr = exp(multi_seas_mod_jnj_fc_df$Forecast+((multi_seas_mod_jnj_fc_df$`Std. Error`)^2)/2)
dummy_1st_fmean_ts <- ts(
    c(c(da_ts[eotr]), as.numeric(multi_seas_e_qer_arr)), 
    frequency = freq, start = end(tr_da_ts)
)
# Label 2: NULL
lines(dummy_1st_fmean_ts)
# Label 3: Forecast Mean
points(dummy_1st_fmean_ts, pch = 1)
# Plot confidence interval (95%)
t_stats = qnorm(p = 0.975, mean = 0, sd = 1)
multi_seas_mod_jnj_fc_res$lower[,2]
dummy_1st_flower_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), multi_seas_mod_jnj_fc_res$lower[,2])), # lower CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
dummy_1st_fupper_ts <- ts(
    exp(c(c(lg_da_ts[eotr]), multi_seas_mod_jnj_fc_res$upper[,2])), # upper CI for log-normal
    frequency = freq, start = end(tr_da_ts)
)
# Label 4: Forecast 95% Lower Bound 
lines(dummy_1st_flower_ts, lty = 2)
# Label 5: Forecast 95% Upper Bound
lines(dummy_1st_fupper_ts, lty = 2)
# Plot original data
orig_plot_ts <- ts(
    da_ts[(eotr - npts + 1):length(da_ts)], 
    frequency = freq, start = time(da_ts)[eotr] - (npts - 1) / freq
)
# Label 6: NULL
lines(orig_plot_ts)
# Label 7: Actual Earning
points(orig_plot_ts, pch = 19)
legend(
    "topleft", 
    legend = c(
        "Actual Earning Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
        "Forecast 95% Upper Bound", NULL, "Actual Earning"
    ), 
#     col = c("black", "red", "blue"), 
    lty = c(1, NA, 2, 2, NA),
    pch = c(NA, 1, NA, NA, 19)
)

1977 ; 1981

CRSP Decile 1 Index¶

da = read.table("../AFTS_data/Ch02/m-deciles08.txt", header = T)
da[1:5,]

crsp_decile_1_rtn = da$CAP1RET
crsp_decile_1_rtn_ts = ts(crsp_decile_1_rtn, frequency = 12, start = c(1970, 1))

jan = rep(c(1,rep(0,11)),39) # Create Jan dummy
m1 = lm(crsp_decile_1_rtn~jan)
summary(m1)

Call:
lm(formula = crsp_decile_1_rtn ~ jan)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.30861 -0.03475 -0.00176  0.03254  0.40671 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.002864   0.003333   0.859    0.391    
jan         0.125251   0.011546  10.848   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06904 on 466 degrees of freedom
Multiple R-squared:  0.2016,	Adjusted R-squared:  0.1999 
F-statistic: 117.7 on 1 and 466 DF,  p-value: < 2.2e-16

plot_time_fig(crsp_decile_1_rtn_ts, main = "Fig. 2.16(a)", xlab = "Year", ylab = "Simple Return")
plot_acf(da = crsp_decile_1_rtn, lag.max = 40, main = "Fig. 2.16(b)")
crsp_decile_1_adj_rtn = ts(m1$residuals, frequency = 12, start = c(1970, 1))
plot_time_fig(crsp_decile_1_adj_rtn, main = "Fig. 2.16(c)", xlab = "Year", ylab = "Adjusted Return")
plot_acf(da = m1$residuals, lag.max = 40, main = "Fig. 2.16(d)")

m2 = arima(crsp_decile_1_rtn, order=c(1,0,0), seasonal=list(order=c(1,0,1), period=12))
summary(m2)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = crsp_decile_1_rtn, order = c(1, 0, 0), seasonal = list(order = c(1, 
    0, 1), period = 12))

Coefficients:
         ar1    sar1     sma1  intercept
      0.1769  0.9882  -0.9144     0.0118
s.e.  0.0456  0.0093   0.0335     0.0129

sigma^2 estimated as 0.004717:  log likelihood = 584.07,  aic = -1160.14

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

par(bg = 'white')
tsdiag(m2,gof=36)

m3 = arima(crsp_decile_1_rtn, order=c(1,0,0), seasonal=list(order=c(1,0,1), period=12), include.mean = F)
summary(m3)

Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”

Call:
arima(x = crsp_decile_1_rtn, order = c(1, 0, 0), seasonal = list(order = c(1, 
    0, 1), period = 12), include.mean = F)

Coefficients:
         ar1    sar1     sma1
      0.1787  0.9886  -0.9127
s.e.  0.0456  0.0089   0.0335

sigma^2 estimated as 0.00472:  log likelihood = 583.68,  aic = -1161.36

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

Sec. 2.9¶

Weekly 1yr and 3yr interest rate¶

da1 = read.table("../AFTS_data/Ch02/w-gs1yr.txt", header = T)
da3 = read.table("../AFTS_data/Ch02/w-gs3yr.txt", header = T)
da1[1:5,]; da3[1:5,]

Exploration¶

r1 = da1$rate
r3 = da3$rate
freq = 52
r1_ts = ts(r1, frequency = freq, start = c(1962, 1))
r3_ts = ts(r3, frequency = freq, start = c(1962, 1))
c1 = diff(r1)
c3 = diff(r3)
c1_ts = ts(c1, frequency = freq, start = c(1962, 2))
c3_ts = ts(c3, frequency = freq, start = c(1962, 2))

par(bg = 'white')
plot(r1_ts, main = "Weekly 1yr or 3yr U.S. Interest Rate (Fig. 2.17)", xlab = "Year", ylab = "Percent")
lines(r3_ts, lty = 2)
legend(
    "topleft",
    legend = c("1yr rate", "3yr rate"),
    lty = c(1, 2)
)

plot(c1_ts, main = "Weekly 1yr U.S. Interest Rate Changes (Fig. 2.20(a))", xlab = "Year", ylab = "Percent")
plot(c3_ts, main = "Weekly 3yr U.S. Interest Rate Changes (Fig. 2.20(b))", xlab = "Year", ylab = "Percent")

par(bg = "white")
plot(x = r1, y = r3, main = "Fig. 2.18(a)", xlab = "1yr rate", ylab = "3yr rate")
plot(x = c1, y = c3, main = "Fig. 2.18(b)", xlab = "diff(1yr rate)", ylab = "diff(3yr rate)")

Fit models¶

Strong serial correlation in ACF of residuals, indicating that there is a unit-root process.
Also conducted ADF unit-root tests. The results show that null-hypothesis cannot be rejected.
- Null-hypothesis is that the time series has unit-root.

m1 = lm(r3~r1)
summary(m1)

Call:
lm(formula = r3 ~ r1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.82319 -0.37691 -0.01462  0.38661  1.35679 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.83214    0.02417   34.43   <2e-16 ***
r1           0.92955    0.00357  260.40   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5228 on 2465 degrees of freedom
Multiple R-squared:  0.9649,	Adjusted R-squared:  0.9649 
F-statistic: 6.781e+04 on 1 and 2465 DF,  p-value: < 2.2e-16

par(bg = 'white')
plot(m1$residuals, type = 'l', main = "Simple Linear Model (Fig. 2.19(a))", xlab = "Year")
plot_acf(da = m1$residuals, lag.max = 36, main = "Residuals ACF Plot (Fig. 2.19(b))")
pacf(x = m1$residuals, lag.max = 36, main = "Residuals PACF Plot")

help(adfTest)

# "lags=12" is chosen from PACF plot
adf_test_res_1 = adfTest(r1, lags = 12, type = c("ct"))
adf_test_res_3 = adfTest(r3, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_1; adf_test_res_3

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -2.4804
  P VALUE:
    0.3749 

Description:
 Tue Feb 27 04:13:44 2024 by user:

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -2.0718
  P VALUE:
    0.5479 

Description:
 Tue Feb 27 04:13:44 2024 by user:

box_test_1 = Box.test(m1$residuals, lag = 12, type = 'Ljung')
box_test_1

	Box-Ljung test

data:  m1$residuals
X-squared = 24787, df = 12, p-value < 2.2e-16

The serial correlation is largely reduced, but there are still some serial correlation in ACF of residuals, indicating that the model is not adequate. We need to model some serial correlation in residuals. The ACF plot shows strong serial correlation at lag-1, so that we use MA(1) to model the error.

m2 = lm(c3~-1+c1) # No intercept
summary(m2)

Call:
lm(formula = c3 ~ -1 + c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42469 -0.03589 -0.00127  0.03456  0.48911 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
c1 0.791935   0.007337   107.9   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06896 on 2465 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.165e+04 on 1 and 2465 DF,  p-value: < 2.2e-16

plot_acf(da = m2$residuals, lag.max = 36, main = "Fig. 2.21(b)")

# "lags=12" is chosen from PACF plot
adf_test_res_c1 = adfTest(c1, lags = 12, type = c("ct"))
adf_test_res_c3 = adfTest(c3, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_c1; adf_test_res_c3

Warning message in adfTest(c1, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”
Warning message in adfTest(c3, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -13.6578
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:11:17 2024 by user:

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -13.933
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:11:17 2024 by user:

box_test_2 = Box.test(m2$residuals, lag = 12, type = 'Ljung')
box_test_2

	Box-Ljung test

data:  m2$residuals
X-squared = 141.84, df = 12, p-value < 2.2e-16

The serial correlations in ACF plots no longer exists. But the Box-Ljung test is still significant. Maybe we can improve the model more to reduce some serial correlation in the residuals.

m3 = arima(c3, order = c(0,0,1), xreg = c1, include.mean = F)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(m3$residuals^2)/sum(c3^2)
summary(m3); rsq

Call:
arima(x = c3, order = c(0, 0, 1), xreg = c1, include.mean = F)

Coefficients:
         ma1      c1
      0.1823  0.7936
s.e.  0.0196  0.0075

sigma^2 estimated as 0.0046:  log likelihood = 3136.62,  aic = -6267.23

Training set error measures:
                        ME       RMSE        MAE MPE MAPE      MASE
Training set -7.953792e-05 0.06782056 0.04843167 NaN  Inf 0.3766059
                     ACF1
Training set 0.0008170501

plot_acf(da = m3$residuals, lag.max = 36, main = "ACF of residuals after taking difference\nand model error as a MA(1) process")

log(length(c1))

box_test_3 = Box.test(m3$residuals, lag = 12, type = 'Ljung')
box_test_3

	Box-Ljung test

data:  m3$residuals
X-squared = 54.422, df = 12, p-value = 2.297e-07

Sec. 2.10¶

Example of consistent estimation for covariance matrix (1yr and 3yr rates)¶

For heteroskedasticity-consistent or heteroskedasticity&autocorrelation-consistent estimators.

da1 = read.table("../AFTS_data/Ch02/w-gs1yr.txt", header = T)
da3 = read.table("../AFTS_data/Ch02/w-gs3yr.txt", header = T)
da1[1:5,]; da3[1:5,]

r1 = da1$rate
r3 = da3$rate
freq = 52
r1_ts = ts(r1, frequency = freq, start = c(1962, 1))
r3_ts = ts(r3, frequency = freq, start = c(1962, 1))
c1 = diff(r1)
c3 = diff(r3)
c1_ts = ts(c1, frequency = freq, start = c(1962, 2))
c3_ts = ts(c3, frequency = freq, start = c(1962, 2))

m1 <- lm(c3~c1)
summary(m1)

Call:
lm(formula = c3 ~ c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42459 -0.03578 -0.00117  0.03467  0.48921 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0001051  0.0013890  -0.076     0.94    
c1           0.7919323  0.0073391 107.906   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06897 on 2464 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.164e+04 on 1 and 2464 DF,  p-value: < 2.2e-16

HC estimator
- GPT: R use white correction in OLS summary

help(vcovHC)

# Applying White's correction
# `type` parameter 'arg' should be one of “HC3”, “const”, “HC”, “HC0”, “HC1”, “HC2”, “HC4”, “HC4m”, “HC5”
hc1_robust_se <- vcovHC(m1, type = "HC1") # HC1 is often recommended; it's a variant of White's correction
hc1_robust_se

help(coeftest)

# Getting the summary with corrected standard errors
summary_with_hc1_robust_se <- coeftest(m1, hc1_robust_se)
coeftest(m1, vcovHC);
coeftest(m1, vcovHC, type = "HC0");
print(summary_with_hc1_robust_se);
dwtest(m1);
jarque.bera.test(m1$residuals);
Box.test(m1$residuals, lag = 33, type = 'Ljung')

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00139139 -0.0756   0.9398    
c1           0.79193231  0.01659956 47.7080   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00138797 -0.0758   0.9396    
c1           0.79193231  0.01634193 48.4601   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00138853 -0.0757   0.9396    
c1           0.79193231  0.01634856 48.4405   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

	Durbin-Watson test

data:  m1
DW = 1.6456, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

	Jarque Bera Test

data:  m1$residuals
X-squared = 1644.6, df = 2, p-value < 2.2e-16

	Box-Ljung test

data:  m1$residuals
X-squared = 230.05, df = 33, p-value < 2.2e-16

HAC estimator
- GPT: Use heteroscedasticity and autocorrelation consistent (HAC) estimator in R

help(vcovHAC)

# Compute HAC standard errors
hac_se <- vcovHAC(m1)
hac_se

# Summary of the model with HAC standard errors
# Doesn't have effect
summary(m1, robust = hac_se)

Call:
lm(formula = c3 ~ c1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42459 -0.03578 -0.00117  0.03467  0.48921 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0001051  0.0013890  -0.076     0.94    
c1           0.7919323  0.0073391 107.906   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06897 on 2464 degrees of freedom
Multiple R-squared:  0.8253,	Adjusted R-squared:  0.8253 
F-statistic: 1.164e+04 on 1 and 2464 DF,  p-value: < 2.2e-16

# Alternatively, you can directly use coeftest() from lmtest to test coefficients with HAC standard errors
coeftest(m1, vcov. = hac_se)

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00158285 -0.0664    0.947    
c1           0.79193231  0.01746723 45.3382   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The following function gives closer answer at the pp100 in the textbook.

coeftest(m1, NeweyWest(m1, lag = 12, prewhite = F)) # From the solution to 2-7

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -0.00010515  0.00149696 -0.0702    0.944    
c1           0.79193231  0.02031114 38.9901   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Including lag in linear regression
- GPT: R linear model using lag operator

help(merge)

# Create a zoo object for the time series
c3_zoo <- zoo(c3)
c1_zoo <- zoo(c1)

# Lag the series by 1 period
c3_lag1_zoo <- lag(c3_zoo, -1, na.pad=TRUE)
c1_lag1_zoo <- lag(c1_zoo, -1, na.pad=TRUE)

# Combine the original and lagged series
c3_c1_zoo <- merge(c3_zoo, c1_zoo, c3_lag1_zoo, c1_lag1_zoo, all = TRUE)
dim(c3_c1_zoo); c3_c1_zoo[1:5,]

  c3_zoo c1_zoo c3_lag1_zoo c1_lag1_zoo
1   0.05   0.08          NA          NA
2   0.05  -0.03        0.05        0.08
3  -0.03  -0.03        0.05       -0.03
4   0.03   0.03       -0.03       -0.03
5  -0.04   0.00        0.03        0.03

help(lm)

# Fit a linear model
model_zoo <- lm(c3_zoo ~ c1_zoo+c3_lag1_zoo+c1_lag1_zoo, data = c3_c1_zoo, na.action = na.exclude)
# View the summary of the model
summary(model_zoo)
dwtest(model_zoo);
jarque.bera.test(model_zoo$residuals);
Box.test(model_zoo$residuals, lag = 33, type = 'Ljung')

Call:
lm(formula = c3_zoo ~ c1_zoo + c3_lag1_zoo + c1_lag1_zoo, data = c3_c1_zoo, 
    na.action = na.exclude)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.44812 -0.03551 -0.00075  0.03408  0.45821 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.686e-05  1.367e-03  -0.064    0.949    
c1_zoo       7.971e-01  7.692e-03 103.632   <2e-16 ***
c3_lag1_zoo  1.766e-01  1.983e-02   8.906   <2e-16 ***
c1_lag1_zoo -1.580e-01  1.745e-02  -9.058   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06785 on 2461 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.8312,	Adjusted R-squared:  0.831 
F-statistic:  4039 on 3 and 2461 DF,  p-value: < 2.2e-16

	Durbin-Watson test

data:  model_zoo
DW = 1.9865, p-value = 0.3649
alternative hypothesis: true autocorrelation is greater than 0

	Jarque Bera Test

data:  model_zoo$residuals
X-squared = 1620.5, df = 2, p-value < 2.2e-16

	Box-Ljung test

data:  model_zoo$residuals
X-squared = 131.6, df = 33, p-value = 1.009e-13

Prove that the equation on pp101 is the $j^{th}$ diagnoal element of the equation (2.49).¶

From equation (2.49) $$ \begin{eqnarray} Cov(\hat{\boldsymbol{\beta}})_{HC} = \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \left[\sum_{t=1}^T \hat{e}_t^2\mathbf{x}_t\mathbf{x}_t'\right] \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \end{eqnarray} $$
The equation on pp101 $$ Var(\hat{\beta}_j)_{HC} = \frac{\sum_{t=1}^T\hat{e}_t^2\hat{v}_t^2}{\left(\sum_{t=1}^T\hat{v}_t^2\right)^2} $$

where $\hat{v}_t$ is the least-squares residual the auxiliary regression $x_{jt}=\mathbf{x}_{-j,t}'\boldsymbol{\gamma}+v_t, t=1,\cdots,T$.

Basically, we want to prove that $$ \begin{eqnarray} \mathbf{e}_j'\left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1} \left[\sum_{t=1}^T \hat{e}_t^2\mathbf{x}_t\mathbf{x}_t'\right] \left[\sum_{t=1}^T \mathbf{x}_t\mathbf{x}_t'\right]^{-1}\mathbf{e}_j = \frac{\sum_{t=1}^T\hat{e}_t^2\hat{v}_t^2}{\left(\sum_{t=1}^T\hat{v}_t^2\right)^2} \end{eqnarray} $$

where $\mathbf{e}_j$ is the unit vector with $j^{th}$ element being $1$ and other element being $0$.

Before proving the general case, we can first check if this is true when all the columns of $\mathbf{X}$ are orthogonal.

We can write $\left[\sum_{t=1}^T\mathbf{x}_t\mathbf{x}_t'\right]=\mathbf{X}'\mathbf{X}$ as a block matrix and use the formula of inversion of a block matrix to prove this.

Sec. 2.11¶

Daily simple return of value-&equal-weighted index from CRSP¶

da_daily = read.table("../AFTS_data/Ch02/d-ibm3dx7008.txt", header = T)
vwretd = da_daily$vwretd
ewretd = da_daily$ewretd

par(bg = "white")
plot(acf(abs(vwretd), lag.max = 300, plot = F), main = "ACF of abs(vwretd) (Fig. 2.22(a))")
plot(acf(abs(ewretd), lag.max = 300, plot = F), main = "ACF of abs(ewretd) (Fig. 2.22(b))")

	date	ibmrtn	vwrtn	ewrtn	sprtn
	<int>	<dbl>	<dbl>	<dbl>	<dbl>
1	19260130	-0.010381	0.000724	0.023174	0.022472
2	19260227	-0.024476	-0.033374	-0.053510	-0.043956
3	19260331	-0.115591	-0.064341	-0.096824	-0.059113
4	19260430	0.089783	0.038358	0.032946	0.022688
5	19260528	0.036932	0.012172	0.001035	0.007679

	date	ibmrtn	vwrtn	ewrtn	sprtn
	<int>	<dbl>	<dbl>	<dbl>	<dbl>
1	19260130	-0.010381	0.000724	0.023174	0.022472
2	19260227	-0.024476	-0.033374	-0.053510	-0.043956
3	19260331	-0.115591	-0.064341	-0.096824	-0.059113
4	19260430	0.089783	0.038358	0.032946	0.022688
5	19260528	0.036932	0.012172	0.001035	0.007679

	date	ibmrtn	vwrtn	ewrtn	sprtn
	<int>	<dbl>	<dbl>	<dbl>	<dbl>
1	19260130	-0.010381	0.000724	0.023174	0.022472
2	19260227	-0.024476	-0.033374	-0.053510	-0.043956
3	19260331	-0.115591	-0.064341	-0.096824	-0.059113
4	19260430	0.089783	0.038358	0.032946	0.022688
5	19260528	0.036932	0.012172	0.001035	0.007679

	date	ibmrtn	vwrtn	ewrtn	sprtn
	<int>	<dbl>	<dbl>	<dbl>	<dbl>
1	19260130	-0.010381	0.000724	0.023174	0.022472
2	19260227	-0.024476	-0.033374	-0.053510	-0.043956
3	19260331	-0.115591	-0.064341	-0.096824	-0.059113
4	19260430	0.089783	0.038358	0.032946	0.022688
5	19260528	0.036932	0.012172	0.001035	0.007679

	date	rtn
	<int>	<dbl>
1	19460228	-0.077922
2	19460330	0.018592
3	19460430	-0.100000
4	19460531	0.209877
5	19460628	0.005128

	Forecast	Std. Error	Actual
Jan 2008	0.007455686	0.05328902	-0.062305
Feb 2008	0.015952488	0.05357340	-0.022040
Mar 2008	0.011701663	0.05357574	-0.010470
Apr 2008	0.009806653	0.05390675	0.051144
May 2008	0.008770886	0.05392026	0.023825
Jun 2008	0.009164644	0.05392029	-0.078625
Jul 2008	0.009432520	0.05392456	-0.013150
Aug 2008	0.009565095	0.05392493	0.011042
Sep 2008	0.009530555	0.05392493	-0.098060
Oct 2008	0.009495151	0.05392499	-0.184726
Nov 2008	0.009477748	0.05392500	-0.085206
Dec 2008	0.009480420	0.05392500	0.021482

	Forecast	Std. Error	Actual
Mar 2008	0.004284404	0.07117457	-0.025971
Apr 2008	0.013560426	0.07237450	0.031213
May 2008	0.015025565	0.07237450	0.032178
Jun 2008	0.014453940	0.07288140	-0.087065
Jul 2008	0.012047087	0.07288140	-0.000994
Aug 2008	0.001806804	0.07288140	0.014066
Sep 2008	0.012211891	0.07288140	-0.120888
Oct 2008	0.005515984	0.07288140	-0.206004
Nov 2008	0.008514008	0.07288140	-0.136597
Dec 2008	0.012792573	0.07339524	0.043083

	0	1	2	3	4	5	6	7	8	9	10	11	12
	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>
1	-0.056	-0.03796	-0.0822	-0.00465	0.017744	0.08207	0.00800	0.01267	-0.030143	-0.07777	0.04877	0.09091	-0.0109
2	-0.474	0.00958	-0.0736	-0.02090	0.001963	0.07722	-0.02876	0.00259	-0.006832	-0.06937	0.03722	0.09383	-0.0242
3	-0.383	-0.34759	-0.0737	0.01595	-0.005528	0.07720	0.02694	0.01197	0.000449	-0.02678	0.02209	0.04284	0.0418
4	-0.177	0.13811	0.3838	-0.02240	0.002324	0.04192	-0.02319	0.01542	-0.004405	-0.02537	0.01849	0.01002	0.0433
5	0.421	0.02873	0.4542	-0.00787	0.000711	0.00254	-0.01402	0.03052	0.011587	0.00422	0.01910	-0.00434	0.0133
6	-0.114	0.21353	0.4490	0.00958	0.202415	-0.00627	-0.00379	0.04032	-0.012936	-0.01233	0.03150	0.01172	0.0277
7	-0.208	-0.25037	0.2431	0.31109	0.167446	-0.03875	-0.00338	0.04292	-0.010090	-0.02600	0.00782	0.01056	0.0373

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	o	o	x	o	o	x	o	o	o	x	o	x	o
1	x	o	x	o	o	x	o	o	o	o	o	x	o
2	x	x	x	o	o	x	o	o	o	o	o	o	o
3	x	x	x	o	o	o	o	o	o	o	o	o	o
4	x	o	x	o	o	o	o	o	o	o	o	o	o
5	x	x	x	o	x	o	o	o	o	o	o	o	o
6	x	x	x	x	x	o	o	o	o	o	o	o	o

	X..newX..i
	<dbl>
nobs	755.000000
NAs	0.000000
Minimum	-0.326128
Maximum	0.229484
1. Quartile	-0.029969
3. Quartile	0.050031
Mean	0.010299
Median	0.008830
Sum	7.776056
SE Mean	0.002319
LCL Mean	0.005747
UCL Mean	0.014852
Variance	0.004060
Stdev	0.063719
Skewness	-0.076466
Kurtosis	1.250917

	X..newX..i
	<dbl>
nobs	248.000000
NAs	0.000000
Minimum	5.468904
Maximum	9.575872
1. Quartile	6.377224
3. Quartile	8.799009
Mean	7.612532
Median	7.664212
Sum	1887.907986
SE Mean	0.081566
LCL Mean	7.451879
UCL Mean	7.773185
Variance	1.649936
Stdev	1.284498
Skewness	-0.063967
Kurtosis	-1.433830

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	x	x	o	o	o	x	o	o	o	o	o	x	o
1	x	x	o	o	o	o	o	o	o	o	o	x	o
2	x	o	o	o	o	o	o	o	o	o	o	x	o
3	x	x	o	o	o	o	o	o	o	o	o	x	o
4	x	o	o	x	o	o	o	o	o	o	o	x	o
5	o	x	x	x	x	x	o	o	o	o	o	o	o
6	x	x	x	x	x	x	o	o	o	o	o	o	o

	X..newX..i
	<dbl>
nobs	14661.000000
NAs	0.000000
Minimum	-0.228997
Maximum	0.087089
1. Quartile	-0.004067
3. Quartile	0.004845
Mean	0.000299
Median	0.000443
Sum	4.382049
SE Mean	0.000074
LCL Mean	0.000153
UCL Mean	0.000445
Variance	0.000081
Stdev	0.009011
Skewness	-1.276328
Kurtosis	33.590044

	Qtr1	Qtr2	Qtr3	Qtr4
1960	-0.34249031	-0.46203546	-0.16251893	-0.82098055
1961	-0.49429632	-0.37106368	-0.08338161	-0.59783700
1962	-0.32850407	-0.26136476	-0.08338161	-0.51082562
1963	-0.18632958	-0.22314355	0.00000000	-0.26136476
1964	-0.08338161	0.00000000	0.21511138	0.00000000
1965	0.14842001	0.26236426	0.37156356	0.22314355
1966	0.23111172	0.32208350	0.62057649	0.44468582
1967	0.42526774	0.46373402	0.60431597	0.62057649
1968	0.42526774	0.72754861	0.85015093	0.81093022
1969	0.77010822	0.88789126	0.99325177	0.81093022
1970	1.02604160	1.22964055	1.30562646	1.28093385
1971	1.28093385	1.46325540	1.46325540	1.39871688
1972	1.58103844	1.61740608	1.61740608	1.48387469
1973	1.71918878	1.76644166	1.88251383	1.66959184
1974	1.79674701	1.85473427	1.93585981	1.76644166
1975	1.93585981	2.04640169	2.05796251	1.81156210
1976	2.04640169	2.18717424	2.11384297	1.92278773
1977	2.25549349	2.32825284	2.25549349	2.16676526
1978	2.47485631	2.48989419	2.49732917	2.18717424
1979	2.64191040	2.56186769	2.69799987	2.30158459
1980	2.78501124	2.68580459	2.77383794	2.45186680

	Qtr1	Qtr2	Qtr3	Qtr4
1979	0.09037218	0.02639844	-0.02343313	0.06359381
1980	-0.03378421	0.00000000	0.00000000	0.00000000

A data.frame: 5 × 4
	year	mon	day	gdp
	<int>	<int>	<int>	<dbl>
1	1947	1	1	237.2
2	1947	4	1	240.5
3	1947	7	1	244.6
4	1947	10	1	254.4
5	1948	1	1	260.4

A data.frame: 10 × 9
	year	mon	day	open	high	low	close	volume	adjclose
	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	1950	1	3	16.66	16.66	16.66	16.66	1260000	16.66
2	1950	1	4	16.85	16.85	16.85	16.85	1890000	16.85
3	1950	1	5	16.93	16.93	16.93	16.93	2550000	16.93
4	1950	1	6	16.98	16.98	16.98	16.98	2010000	16.98
5	1950	1	9	17.08	17.08	17.08	17.08	2520000	17.08
6	1950	1	10	17.03	17.03	17.03	17.03	2160000	17.03
7	1950	1	11	17.09	17.09	17.09	17.09	2630000	17.09
8	1950	1	12	16.76	16.76	16.76	16.76	2970000	16.76
9	1950	1	13	16.67	16.67	16.67	16.67	3330000	16.67
10	1950	1	16	16.72	16.72	16.72	16.72	1460000	16.72

A data.frame: 7 × 13
	0	1	2	3	4	5	6	7	8	9	10	11	12
	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>	<I<chr>>
1	0.0694	-0.0337	-0.01280	-0.002520	5.63e-04	-0.01801	-0.01383	0.000595	-3.82e-03	0.002069	-0.006495	0.02531	0.000255
2	0.4254	-0.0550	-0.00613	-0.004956	-6.03e-05	-0.01229	-0.01443	-0.001431	-2.37e-03	-0.001021	-0.000153	0.02540	-0.000565
3	-0.1919	-0.0147	0.01549	-0.000335	1.72e-03	-0.01139	-0.01638	-0.000969	5.14e-03	-0.000838	-0.000215	0.02057	0.006086
4	-0.2794	-0.1311	-0.00067	-0.007149	1.36e-02	-0.00567	-0.01635	-0.000868	6.89e-04	-0.001961	-0.000320	0.02061	0.004844
5	0.0655	-0.0129	0.00773	-0.194153	-5.89e-03	-0.00823	-0.01319	-0.002503	2.39e-05	-0.001603	-0.000660	0.01698	0.009937
6	0.0085	-0.1240	-0.34471	-0.230192	2.87e-02	-0.04058	-0.00963	0.003877	-4.63e-03	0.004135	-0.008189	0.01272	-0.007121
7	-0.4457	0.2788	-0.31801	0.262998	-3.33e-01	0.11436	0.00744	-0.008721	-4.15e-03	0.002669	-0.001249	0.00202	-0.001039

A data.frame: 6 × 1
	V1
	<dbl>
1	0.71
2	0.63
3	0.85
4	0.44
5	0.61
6	0.69

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	o	o	x	o	o	x	o	o	o	x	o	x	o
1	x	o	x	o	o	x	o	o	o	o	o	x	o
2	x	x	x	o	o	x	o	o	o	o	o	o	o
3	x	x	x	o	o	o	o	o	o	o	o	o	o
4	x	o	x	o	o	o	o	o	o	o	o	o	o
5	x	x	x	o	x	o	o	o	o	o	o	o	o
6	x	x	x	x	x	o	o	o	o	o	o	o	o

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	x	x	o	o	o	x	o	o	o	o	o	x	o
1	x	x	o	o	o	o	o	o	o	o	o	x	o
2	x	o	o	o	o	o	o	o	o	o	o	x	o
3	x	x	o	o	o	o	o	o	o	o	o	x	o
4	x	o	o	x	o	o	o	o	o	o	o	x	o
5	o	x	x	x	x	x	o	o	o	o	o	o	o
6	x	x	x	x	x	x	o	o	o	o	o	o	o

	Qtr1	Qtr2	Qtr3	Qtr4
1979	0.0916983	0.1096627	0.1096627	0.1096627
1980	0.1142420	0.1161564	0.1161564	0.1161564

	Qtr1	Qtr2	Qtr3	Qtr4
1979	0.14664510	-0.09508059	0.12869720	-0.08626034
1980	0.02869049	-0.01916394	-0.04809882	0.07444413

	Forecast	Std. Error	Actual
1979 Q1	0.09037218	0.0916983	0.14664510
1979 Q2	0.02639844	0.1096627	-0.09508059
1979 Q3	-0.02343313	0.1096627	0.12869720
1979 Q4	0.06359381	0.1096627	-0.08626034
1980 Q1	-0.03378421	0.1142420	0.02869049
1980 Q2	0.00000000	0.1161564	-0.01916394
1980 Q3	0.00000000	0.1161564	-0.04809882
1980 Q4	0.00000000	0.1161564	0.07444413

	Qtr1	Qtr2	Qtr3	Qtr4
1979	2.585637	2.627074	2.611076	2.364514
1980	2.729193	2.770630	2.754632	2.508070

	Qtr1	Qtr2	Qtr3	Qtr4
1979	0.09169838	0.09697557	0.10198005	0.10675017
1980	0.14042389	0.14977390	0.15857356	0.16690994

	Qtr1	Qtr2	Qtr3	Qtr4
1979	2.641910	2.561868	2.698000	2.301585
1980	2.785011	2.685805	2.773838	2.451867

	Qtr1	Qtr2	Qtr3	Qtr4
1979	2.405912	2.437005	2.411198	2.155288
1980	2.453968	2.477078	2.443833	2.180933

	(Intercept)	c1
(Intercept)	1.928016e-06	-4.212675e-07
c1	-4.212675e-07	2.672754e-04

Table of Contents

Sec. 2.2¶

Box-Ljung test (IBM and CRSP Value-Weighted Index)¶

Find IBM monthly and yearly expected return¶

Sec. 2.4¶

AR(3)¶

AIC¶

AR(3) example¶

Parameter estimation for value-weighted index¶

Calculate the actual yearly return¶

Model checking¶

AR(3) forecasting example¶

Find Volume-weighted monthly and yearly expected return¶

Sec. 2.5¶

MA(9) fitting with fixed coefs example (Sec. 2.5.3)¶

MA(9) forecasting (Sec. 2.5.4)¶

Find Equal-weighted monthly and yearly expected return¶

Sec. 2.6¶

3M log price (Sec. 2.6.3)¶

Find 3M monthly and yearly expected return¶

Sec. 2.7¶

3M log price¶

Unit-root test (U.S. quarterly GDP)¶

Full ARIMA workflow (S&P 500)¶

Unit-root testing¶

Order determination and parameter estimations¶

Box-Ljung test¶

Plot forcasting¶

Find S&P500 monthly and yearly expected return¶

Sec. 2.8¶

J&J quarterly earnings¶

Plot the acf of JnJ quarterly earnings¶

Apply airline model to the JnJ data.¶

Box-Ljung test¶

Forecast¶

Forecast: Directly apply multiplicative seasonal model¶

CRSP Decile 1 Index¶

Sec. 2.9¶

Weekly 1yr and 3yr interest rate¶

Exploration¶

Fit models¶

Sec. 2.10¶

Example of consistent estimation for covariance matrix (1yr and 3yr rates)¶

Prove that the equation on pp101 is the $j^{th}$ diagnoal element of the equation (2.49).¶

Sec. 2.11¶

Daily simple return of value-&equal-weighted index from CRSP¶

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	o	o	x	o	o	x	o	o	o	x	o	x	o
1	x	o	x	o	o	x	o	o	o	o	o	x	o
2	x	x	x	o	o	x	o	o	o	o	o	o	o
3	x	x	x	o	o	o	o	o	o	o	o	o	o
4	x	o	x	o	o	o	o	o	o	o	o	o	o
5	x	x	x	o	x	o	o	o	o	o	o	o	o
6	x	x	x	x	x	o	o	o	o	o	o	o	o

	0	1	2	3	4	5	6	7	8	9	10	11	12
0	x	x	o	o	o	x	o	o	o	o	o	x	o
1	x	x	o	o	o	o	o	o	o	o	o	x	o
2	x	o	o	o	o	o	o	o	o	o	o	x	o
3	x	x	o	o	o	o	o	o	o	o	o	x	o
4	x	o	o	x	o	o	o	o	o	o	o	x	o
5	o	x	x	x	x	x	o	o	o	o	o	o	o
6	x	x	x	x	x	x	o	o	o	o	o	o	o

	(Intercept)	c1
(Intercept)	2.505407e-06	-1.306652e-06
c1	-1.306652e-06	3.051041e-04